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ABSTRACT 



In testing the hypothesis that there is no itcrotone 
trend in a gamma renewal process, the use of the statistic 



where 



Y 

J 




/ Y 



1J 






J 

Y = £ X 

1 J i= 1 i 



and 



J 

Y = J S , 
2 J i-1 i 



is investigated.. The mean and variance of Y 

J 



is developed 



as a function of J and it is shown that Y is asymptotically 

J 

normal as J — > co for the gamma renewal process. A 
high-speed, theoretically exact gamma pseudo-random variate 
generator is developed, tested and compared kith other known 
techniques. The . generator is then used to ofctain the 

distribution of Y through digital computer simulation for 

J 

small and moderate values of J. 
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I. INTRODUCTION 



A renewal process is a model for a series of events 
which occur at random times; the times between each two 
successive events are independent and identically 
distributed positive random variables. The time between 

event i-1 and event i is denoted by X , i=1,2,..., while 

i 

the random variable S represents the time to the i-th 

i 

event. When zero is taken as the time of initiation of the 
series of events then 



S 

i 



1 

Z x 

n= 1 n 



(i 1 r 2 /.».)» 



An equivalent representation of a series of events is 
the counting function N(t) which represents the number of 
events occurring in the interval (0,t], The equivalence of 
the two representations is expressed in the fundamental 
ident ity 



N (t) < n if and only if S > t (n = 1 , 2 r . . . ) . 

n 



Thus P{N(t) < n} = P {S > t} , where the notation P{ } 



represent 
Two 
are the 
expected 
the rate 



n 

s the probability of the event within the brackets, 
important characteristics of any series of events 
so-called renewal function M(t) f which is the 
number of events up to time t f and its derivative 
function: 



M (t) = E[ N (t) ], 

X (t) = d M(t) . 

dt 
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A stationary renewal process is one in which N (t) has 
stationary increments (see Cox and Lewis [6], Chapter 4, for 
details). In this case, the renewal function is 

M (t) = t / u, 

and X (t) = 1/u, where u is the mean time let ween events. 

For an arbitrary series of events there are many 
alternatives to the renewal process. Successive inter-event 
times may not be independent or may depend on the entire 
history of the process, or the distr ibutions of the times 
between events or the increments of K (t) may change in some 
manner as the process continues. When these distributional 
changes involve the expected values of the inter-event times 
or the value of X (t) , then a trend is said to exist in the 
process. The trend way be on the serial number of the 
event; for example, a radio set may tend to fail more 
freguently as a function of the total number of failures 
which have occurred. Alternatively, the trend may be 
essentially independent of the number of events and reflect 
basically the elapsed time t; such a trend cn time might be 
observed in the arrival of jobs at a computer facility over 
a 24-hour day where the job arrival rate depends mainly on 
the time of day and not on the number of jobs submitted. A 
trend on time is usually expressed through a non-linear time 
dependence in M(t). 

Trends may be increasing, decreasing, cyclic cr some 
combination of the three. The various statistical 
procedures designed to test whether an arbitrary series of 
events has a ti~end of any type are guite specific as to the 
type of process and the trend model assumed. For a more 
complete discussion of tests for trend in series cf events 
see Cox and lewis [6], Cox [7] and Lewis [18j. 

The simplest form of renewal process is the Foisson 
process in which the inter-event times have the negative 
exponential distribution. A more general process is the 
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ga mma 
gauitna 



renewal process in which the 
random variables with density fa 



X ' s are 
i 

notion 



independent 



f (x) 



rW 



k- 1 k - X x 
X x e 



(k > 0) . 



In this pararaeteriza.tion, X is a scale parameter while k 
characterizes the distribution and is called the shape 
parameter. The function T(k) is the complete gamma 
function or integral 



r<k) 




k-1 — x 
x e dx . 



In testing for trend on serial number in the gamma 
renewal process the use of the statistic 



where 



and 



Y = Y / Y , 
J 2 J 1 J 



Y = f X 
1 J i-1 i 



Y = J S 
2 J i= 1 i 



is investigated in this paper. h test 

developed, and the distribution of 

through sinulation for small values of 
theory results as J — > co . 



for trend using Y is 

iJ 

Y is investigated 
J 

J and through normal 



8 



II. THEORETICAL RESULTS 



A. BACKGROUND 



Cox [5] proposed the use of the statistic 



B = 2 t 

i= 1 i 



as a test for trend in a Poisson process, where t ,t 



, t 



12 J 

are the normalized times to events in a fixed observation 
interval of length T and J is the number of events occurring 
in (0,T], that is, the observed value of N (T) . Bartholomew 
[4] dealt with the statistic 



m = 1 B 
‘J 

and developed results on the asymptotic relative efficiency 
of tests for trend based on this statistic against several 
different trend models. Lewis [18] showed that the test 
based on m is not a consistent test for stationary Poisson 
processes against stationary gamma renewal alternatives, 
although it had been used for that purpose. 

The test based on the statistic B is a conditional test, 
the conditioning being on the observed value J of the random 
variable N (T) . This is because the test is designed to 
detect trend regardless of the rate of occurrence of events. 
The test arises from the non- homogeneo us Poisson process 
with rate function 



X (t) = exp[a + bt} = X e 



bt 
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In this case it is desired to test whether b = 0, that is 
whether there is a trend. It is possible to obtain the 
likelihood function L ( X , b) of the observations for this 
model and from Hey man-Pearson theory conclude that the test 
based on E is the uniformly most powerful conditional test 

H : b=0, X >0 
0 

H : b*0, X >0 
1 

for every value of T. For details, see Cox and Lewis [6], 
chapters 2 and 3. 

For the foregoing model, the parameter X is a nuisance 
parameter since it has no effect on the trend; it merely 
controls the base rate of the process. It is eliminated in 
the B test by conditioning on J = N (T) , this being a 
sufficient statistic for X for any fixed value of b. 

There is an analogous situation when testing for trend 
on serial number and the process is observed for a fixed 
number of intervals J. Here, the trend model is one of 
exponentially distributed intervals with 

bi 

E[ X ] = exp {a + bi} = Xe (i=1,2,...). 

i 



The sufficient statistic for X is now S 

J 



Y so that the 
1 J 



test is now conditioned on the observed value of Y and the 

1 J 

statistic analogous to B is 



Y 

2 J 



J J 

IS, = ? (J + 1-i) X.. 

i=1 l i=1 l 



Asymptotically , 
conditioning on a 
fixed number of ev 
convenient for t 



there is no 
fixed time interval 
ents but the latter 
he gamma renewal 



difference 


bet 


ween 


and conditio 


ning 


on a 


turns out to 


be 


more 


process. In 


this 


case 
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performing a test for trend analogous to the Poisson process 
test based on B requires the determination of the 

conditional distribution of Y given Y 

2J 1J 



B. THE GAMMA RENEWAL PROCESS 

In the gamma renewal process, it is both analytically 
and computationally simpler to test for trend on serial 
number over a fixed number of observations than to test for 
trend in time over a fixed time period. The computational 
advantage of the former test arises from the ease of 
computer simulation of a fixed number of events as compared 
to the difficulties inherent in continuous time simulation 
control. Analytically, the gamma renewal process has 
certain theoretical properties which simplify consideration 
of the serial trend test. 

For the gamma renewal process with shape parameter k and 
the assumed trend model (Cox [7]) 



E[ X ] = exp {a + bi} = X e 
i 



bi 



(i = 1 , 2, . . .) 



it is possible to set up the likelihood function L(k,b, X ) » 
Explicit results, however, cannot be obtained for estimating 
b or testing whether b=0; it is therefore usual to use the 
Poisson test based on B in this case. The test will 
probably be relatively powerful for k near one but not when 
k is small as it uses very little information about the 



individual X 's. In fact, the test is 
i 



based only on the 



centroid of the times to events and not on times between 
events. 

Elimination of the nuisance parameter X cannot be 
accomplished for the gamma renewal process by conditioning 



on S 



as 



it could in the Poisson case, 



since S is not a 
J 



J 

sufficient statistic for X for all values of b. However, 
since X is a pure scale factor the ratio statistic 



Y 

J 



Y / Y 
2J 1 J 



is free cf X for any b and k and i/s thus a reasonable 
alternative statistic on which to base a test which is free 
of the scale factor. 

The distributions of the two test statistics 
(conditional and ratio) are equivalent in the gamma case 
under the null hypothesis (b=0) because of a result which 
characterizes the gamma distribution. Laha [16] and Lukacs 
£21] have established that, for a sequence of independent 

and identically distributed random variables X ,X , ... ,X , 

12 n 

the statistics 



J 

J (J+1-i) X / 
i = 1 i 

and 



J 

5 X 
i=1 i 



J 

ix 

1 = 1 : 



are independent if and only if the 




distribution . 
renewal process 

E {Y < 

J 



Thus, Y is independent of 
J 

. Using this result. 



= P { Y 


/ 


Y 


< z} 




2 J 


1J 




= P (Y 


/ 


Y 


^ z | 


Y 


2 J 




1J 




1J 


= P{Y 


< 


zy 1 


Y 


= y} 


2 J 




Id 



have the gamma 

Y for the gamma 

1J 

y} 
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so that the conditional and ratio statistics have the same 
distribution. Furthermore, the gamma renewal process is the 
only one for which this is the case. 

A further result of the above characterization and the 
probability relationship is that, if the joint distribution 

of Y and Y is asymptotically bivariate normal (a result 
1J 2J 



which is proved in the following Section) , then Y must also 

be asymptotically normal. This follows since the 

conditional distribution of Y given Y will be 

2 J 2 J 

asymptotically normal under these conditions. Thus, 
although the asymptotic distribution of a ratio statistic 



such as Y can be very difficult to obtain in general, an 
J 

expression can be found for the gamma case. 



Trend tests using Y for other types of renewal process 

J 

could and should be considered but the gamma process uas 
chosen deliberately to take advantage of the above 
simplifying distributional results. The gamma process would 
probably not be the simplest case to consider under the 
trend alternative, however. 

In the succeeding sections, results on the asymptotic 



distribution of Y and the mean and variance of Y for small 

J J 

values of J are developed. In addition, a description of a 

computer simulation experiment which was performed to obtain 

the distribution of Y for small values of J is described 

J 

and its results presented. The power of the Y test for the 

J 

gamma process has not been investigated nor has its power 
relative to nonp arametric tests or tests using data 
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transformations been examined. 



C. ASYMPTOTIC DISTRIBUTIONS 



The statistics 
defined by 



Y 

1 J 



a nd 



Y for a renewal process are 
2 J 



Y 

1 J 



J 

Z X 
1=1 i 



e 



J 

Y = ? S 

2J i = 1 i 



J 

.2 (J+1-i) x., 

i=l i 



where X ,X r ... ,X are the process inter-event times and 
1 2 J 

2 

are assumed to have finite mean u and finite variance a 



Obviously, Y and Y are dependent random variables 
1 J 2 J 

with respective means and variances 



ra 

1 



E[ Y ] 
1 J 



E[ .2 X ] 
1=1 1 



Ji E[ V 



J E[ X ] 
J u. 



2 

s = Var[Y ] 

1 2J 

= Var[ i X. ] 

i = 1 l 
J 

= £ Var[ X ] 

i= 1 i 

= J Var[ X ] 
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2 

J a » 



To 

and 



E[Y 2J ] 



E[ 7 (J+ 1-i) X ] 

i = 1 i 



.2, (J+1-i) E[ x. 3 
i = 1 i 



E[ X] • [ 2 (J+1) 
1=1 



J 

il , 1 1 



2 



s = 



Var[Y ] 
2 J 



Var[ 5 (J+1-i) X ] 
i =1 i 

J 2 

> (j + i-i) Var[ X ] 
i = 1 i 

J 2 

Var[ X ] 2 i 

i=1 






find the correlation coefficient for the statistics Y 



1 J 



Y it is necessary to determine the first joint moment: 
2 J 



= E[ Y Y ] 

12 1 J 2J 

= E[ .2 X. * .2 (J+ 1-i) X. ] 

1 = 1 1 1 = i 1 

= E[.i 2 (J+l-k) xx ] 

i=1 k=1 i k 
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J 2 

= E[ 7 (CM- 1 - i) X 
12 i=1 i 

J-1 J 

+ J 7 (2J+2-i-k) x x ] 

i=1 k=i+ 1 i k 

J 2 

= J (J+ 1-i) E[ X. ] 

1=1 l 

J-1 J 

+ 7 | ( 2 J+2-i-k) E[ X X ] 

i = 1 k=i+'l i k 

2 J 2 J-1 J 

= E[X ] .Z i + (E[X]} 7 , Z |i + k) 

i = 1 i=1 k= i+ 1 



After some simplification this becomes 



2 2 2 

m = J(J+1) Ef X ] + J (J -1) u . 
1 2 2 2 — 



Then the covariance is 



u = E [ ( Y - m ) (Y - m ) ] 

11 1J 1 2J 2 



= m -mm 
12 1 2 



= U3 + 11 a 



and tlie correlation coefficient is 



P 



u / s s 
11 12 

2 

JJJ+11 r / 






/3I3±IEI23±3E o 2 



= /5 
2 



+T 
23TT 



To obtain asymptotic results for Y as J — > » , it is 

iJ 

necessary to show that the joint distribution of Y and Y 

1 J 2 J 

is asymptotically a bivariate normal distribution. To do 
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this it is necessary to apply the Central Limit Theorem to 
the two-dimensional vector random variable 



Z = (Y Y ) 
J 1J 2 J 



Wald and Wolfowitz [29] have shown that Z is asymptotically 

iJ 



normal if and only if 



L = 1 Y + 1 Y 
J 1 IJ 2 2 J 



is also asymptotically normal for every choice of 1 and 1 . 

1 2 

To apply the Central Limit Theorem to L , define 

J 



L* = L - E[L 1 
J J J 

= 1 .2 X + i N + i-i) X. - 

1 1=1 l 2 i=1 l 



E[L ] 
J 



z [1 + 1 (J+ 1-i) ] (X - u) 

i= 1 1 2 i 

la, (X. - u) . 

i= 1 l l 



a =1 +1 (J+1-i) . 

i 12 



F (v) = P [a (X - u) < v} 
i i i 



s = Var[ L ] 
1 J 



= V a r[ L * ] 

J 

J 2 

= Var[ X - u ] 2 a 

i i= 1 i 
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o 



2 



2 J 
? a 
i= 1 i 



Then L will be asymptotically normal if and only if the 
J 

Lindeberg condition is satisfied for every e > 0 (Cramer 

[ 8 ]) : 



( 2 - 1 ) 



lira 
J — > 



1 

s* 

1 



J 

ill 



c,f 



2 

v dF (v) } 
v I >es i 

1 



= 0, 



After making the substitution v = a z the integral becomes 

i 



f v dF (v) = f a z dF (a z) . 

‘ I v 1 > es i 0 | a z | >e s i i i 

1 i 1 

This can be simplified by first noting that 



F (a z) = P {a (X - u) < a z} 
i i i i i 



= P [X - u < z) 
i 



= P {X < z+u) 
i 

= F (z+u) , 



where F (x) is the common distribution function of the X *s. 

i 

Secondly, if 



| a z | > e s 
i 1 



then 



| z. | > e 



s i 7 



> e s / 1 1 
1 1 



+ 



1 (J + 1-i) | 
2 
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> e s / ( 1 1 | + 1 1 (J + 1-i) | ) 

1 1 2 

> e s / ( 1 1 | + 1 1 I J) - 

1 1 2 



We can thus assume, without loss of generality, that 1 >0 

1 

and 1 >0. Then, 

2 



2 2 2 
v dF (v) < a f zdF (z + u) . 

|vj>es i i |z|>es /(I +1 J) 



J « 2 

1 2 {J v dF (v)} 

‘s? i=1 M | v | >es i 

^ ^ j 2 2 

< 1 [ y a ] f z dF (z+u) 

i= 1 i J | z i>es /(I +1 J) 

1 112 



( 2 - 2 ) 



1 f 



" | z | >es / (1 +1 J) 
1 1 2 



z dF (z + u) 



Since by hypothesis the X 1 s have a finite second moment and 

i 



since 



lim 
J — > 



es /(I +1 J) 
<» 1 12 



> (e/1 ) lim (s /J) 

2 J — > «> 1 

J 2 1/2 

> (eo /I ) lim 1 { ? a } 

2 J— >09 J l iHi i 

J 2 1/2 

> (eg /I ) lim {.2, 1 [1 +1 (J+1-i) ] } 

2 j — > co i-i 12 



the integral (2-2) approaches zero and so the Lindeberg 
condition (2-1) is satisfied. 
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It has thus been shown that Y and Y are 

1 J 2J 

asymptotically jointly normal. From the known properties of 
the two-dimensional normal distribution (Gnedenko [10]), the 

conditional distribution of Y given Y must also be 

2 J 1J 

asymptotically normal with mean and variance 



E[ Y |Y =y ] = m + (P s / s ) (y - m ) 
2 J 1 J 1 2 2 111 



OJil+JL u + iil+11 (y - J u) 



J + 1 y 

1 



Var[ Y | Y =y ] 
2J 1J 1 



2 2 

V 1 - p » 



intiii 



2J+.1]_ 



2 2 
= J_[J -1]_ c . 
72 •* 



6 J-t-1 ) 

£T 2777 



From the results of the previous Section, 
asymptotically normal with mean and variance 



Y is also 
J 



E[ Y ] = E[ Y 
J 2 J 



/ * 



1 J 




= J + 1 . 

~' 2 ~ 



VarfY^] 



= V a r [ Y / Y | Y = y ] 

2 J 1J 2J 1 

2 2 

= -11 ( ^ / y ) - 

12 1 



The foregoing result then suggests the following 
procedure to test for trend in a gamma renewal process. The 



20 



observed values of Y 



can be used to accept or 



and Y 

1 J 2 J 

reject the null hypothesis of a trend-less process at any 
desired significance level a using the conditional normal 
distribution given above. Rejection of the no-trend 

hypothesis will take place if the observed value of Y is 
greater than the 1 -a /2 quantile of the distribution of Y 



(for increasing trend) or if Y is less than the a/2 

quantile (decreasing trend) . 

For the distribution of the sum of a series of 
independent exponential random variables to become 
approximately normal,. however, requires over 20 terms; the 
gamma renewal process with shape parameter less than one 
reguires several terms just to produce a single exponential 
variable. Thus, on the order of 100 events could be 
reguired to use the asymptotic test. For a small sample 
(roughly, one with fewer than 20/k events) it may not be 
possible to use a normal approximation to the distribution 
of the test statistic. 

If sample size is so small that the asymptotic result 
cannot be used, it is then necessary to develop the 

distribution of Y in order to use this test for trend. If 

J 

the distribution cannot be determined analytically, computer 
simulation must be resorted to; it is still possible, 
however, to use the probability relationship between the 

distribution of Y and Y given Y to obtain the moments 

J 2 J 1 J 



of Y . 
J 
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D. THE WEAN AND VARIANCE OF Y FOR SHALL SAMPLES 

J 



From the foregoing, it can be concluded that in a 
renewal process 



Ef Y ] = E[ Y / Y ] 
J 2J 1 J 



= E[Y ] / E[Y ] 
2 J l J 



since Y /Y is independent of Y . Thus, 
2 J 1J 1J 



E [ Y ] = m / m 
J 2 1 



= J+ 1 . 
“ 2 “ 



Similarly, 



E[Y ] 

iJ 



2 2 

= E[ Y / Y ] 

2 J 1J 

■ b [ 4 ]/ e [ 4 ] 

2 2 2 2 

= ( s 0 + %) / ( s . + ®> 

2 2 11 

2 2 2 2 

d c * + >J2u^ 

2 2 

= 1J+1JL lllltll c + 3 J ( J + 1 ) u . 

T2 o'^"+ J u 2- 



Introducing the mean and variance of the gamma 
variable , 



u = k / \ 



gamma 



random 
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2 2 
c = k / X 



E[ V = 

= J+1 2{2J+1) +_3JlJ±lLk . 

“ T2 ~ kJ + T 

The final result for the variance is then 

2 2 
Var[ Y ] = E[ Y ] - {E[ Y ]} 

U J u 

= J- 1 J+ 1 . 

72“ TT3T1 

Note that for the Poisson process (for v^/hich k = 1.0), 

VarfY ] = J-1 , which is also the result for the sum of J-1 

J T2~ 

uniform random variables- When k^l.O then, approximately. 



Var[ Y ] = J-1 C (X) , 
J T2“ 



2 

where C (X) is the coefficient of variation (which is 1/k 

for the gamma random variable.) Thus, -the greater 
variability in the inter-event times is reflected directly 

in the increased variability of the statistic Y . From the 

J 

point of view of a trend test, the variability of Y is 

J 

increased to allow for the possible confusion of an actual 
trend with a sequence of long intervals in the stationary 
process. As mentioned previously, this indicates that more 
powerful tests for trend could be developed which use more 
information about the interval distribution than just the 
coefficient of variation. 

The succeeding sections discuss determination of the 
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distribution of Y for small samples through digital 

cf 

computer simulation. Values of the predicted mean and 

variance of Y were calculated for each simulation run and 
J 

the results are tabulated as follows (note that the mean 
does not depend on k) : 



J 


MEAN 






VARIANCE 








k = .10 


k = .25 


k = .50 


k = .75 


k = 1. 


10 


5.5 


4.125 


2. 375 


1.375 


. 971 


.750 


30 


15 .5 


18.729 


8. 814 


4.682 


3.188 


2.417 


50 


25.5 


34.708 


15.428 


8.010 


5.409 


4.083 


100 


50.5 


75.750 


32. 048 


16.338 


1 0.964 


8.250 
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III. GE NERATING GAMMA VARIATES FOR S IMUL ATION 

Although there were several known methods of generating 
gamma random variates (see Section. III.E below) , all were 
judged deficient for large-scale computer simulation use; 
some are statistically inadequate for' large samples and all 
require too much computer time. A new method was therefore 
developed to overcome these difficulties. 



A. DESCRIPTION OF ALGORITHM 

Marsaglia's r ectangle-wedge- tail methods for generating 
normal f 14,23] and exponential [22] random variates were 
modified to produce unit gamma variates with shape parameter 
k less than one. The basic approach is to select with 

respective probabilities P ,P , ... ,P one of a series of 

12 n 

random variables with distribution functions F (x),F (x) , 

1 2 

... ,F (x) . If x is the generated random number then 
n 

F(x) = P F (x) + P F (x) + ... + P F (x) . 

11 2 2 n n 

Taking derivatives, 

f (x) = P f (x) + P f (x) + ... + P f (x) . 

112 2 n n 

The method is then equivalent to a geometric decom position 
of the density function of the variate to be generated. For 

efficiency, the P 's and F *s are chosen so that most of the 

i i 

time the generated variate is selected from a simple 
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distribution (preferably a uniform distribution). 

Gamma variates with any desired scale parameter are 
obtained by multiplying a unit gamma variate (that is, one 
with scale parameter one) by the reciprocal of the desired 
scale parameter \ „ Thus, the method developed generates 
only unit gamma variates, which have density function 



f (x) = 



1 

7W 



k-1 -x 

x e 



The unit variates can then be scaled to produce the desired 
result. 

It is difficult to apply the decomposition technique to 
the gamma distribution because the density function has a 
singularity at the origin when the shape parameter is less 



than one. A low value of x (called x ) is thus selected as 

1 

the starting point for the decomposition; values of the 

gamma variate less than x are generated with probability P 

1 0 

by means of a series approximation to the inverse of the 
incomplete gamma function (see Subsection III.B. 1.) A 

series of x 's is then determined by the relations 
i 



h = .25 x 

1 1 

(3^1) h = 2 h 

i i- 1 

x = x + h (i=2, 3, . . . , I) . 

i i-1 i-1 

This procedure continues until the area under the density 

function curve from zero to x is greater than 0.999; the 

value N is then defined to be 1-1. Dependent upon k, the 
parameter N varies between 20 and 100. Values of the 
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variate greater than x are generated by an iterative 

N*1 

procedure (Subsection III. B* 2) with probability p 

N + 1 

The density function has thus been decomposed into N+2 

regions; a "head" region below x f a "tail" region above 

1 

x and a series of N vertical strips. Each strip is 

N + 1 

divided into a large rectangle and a small wedge, as shown 
in Figure 1. The probability for the strip as a whole can 




Figure 1. Decomposition of strip i. 

be found by using the incomplete gamma function, 

« x k-1 -y 

s<kiI) = rra- J 0 y e dy ' 

which is readily evaluated by an infinite series [ 1 ] 

co n n + k 

g (k; x) = 1 5 j[- 1} x 

tW n=0 nT^ii+Tc - } 

suitably truncated for computational purposes. 

If P is the probability for the strip as a whole, p 
i i 



the probability of the rectangle and g the probability of 

i 
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the wedge then the following relations exist: 



P = g (k;x . ) - g (k; x . ) , 

i i+I 1 



p = h . f (x . ) , 

i l i+I 



q = P - p . 
i i i 



For small values of i, g » p because of the extreme 

i i 

steepness of f (x) near the origin; this accounts for the 

increasing values of h in (3-1) . For i greater than 5, 

i 

this relationship reverses so that most of the area under 
the f (x) curve lies within the rectangles. In fact, it was 
found empirically that 

N 

£ p = 0.72 
i- 1 i 

for most values of k in the range 0.1 to 1.0. Thus, the 
gamma random variable can be sampled from a uniform 
distribution over 70 per cent of the time. Presumably, an 
increase in this probability would result from choosing each 

h optimally, but. the gains to be made are small, 
i 



1 . Bina ry Se ar ch Scheme 



Kith the decomposition complete, there remain two 
problems: efficient selection of one of some 100 sampling 

distributions and methods for sampling from each of the four 
different types of distribution (head, tail, rectangle and 
wedge) . 

To formalize the selection problem, suppose that 

events f ,f , ... ,f are to be selected with probabilities 

1 2 M 
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p ,p , ... ,p , respectively. In the gamma generator, M = 

1 2 M 

2K+2. A cumulative probability vector P = (P ,P , ... ,P ) 

12 M 

is defined by 



p . = Z.P. 

x g = 1 3 



(i=1,2,.. .,M) . 



If a random number U uniformly distributed on (0,1) is 
generated, the selection of an event can be carried out by 



comparing U serially v/ith P ,P , ... and stopping when, for 

1 2 



the first time, U < P for some n; event n is then selected. 

n 

This method requires on the order of M/2 comparisons and 
would be far too slow to solve the selection problem. 

Substantial saving of comparisons results from 
applying the binary search method to the problem of finding 

n such that P < U < P . For this scheme, an interval of 

n- 1 n 

uncertainty is defined by selection of the indices i and j 



so that P SUSP at all times; initially, i is taken to 

i j 

be zero (with P =0) and j to be M. U is then compared 

0 

with P and the interval of uncertainty adjusted 

[ (i + j)/2] 

according to the result. This procedure is continued until 
the difference i-j is equal to 1; then n = j. It is not 



difficult to show that log M comparisons are neede 

2 

n, so that this method is considerably faster 
sequential search described previously. 

A further saving of comparisons is still 

however. Since the exact distribution of the p 's 

i 



d to find 
than the 

possible , 
is known. 
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the search procedure can be "skewed” so that more probable 
events can be found in fewer comparisons, thus minimizing 
the average number of comparisons. This can be done by 

comparing U with the value of P which is closest to the 

k 



average of P and P and proceeding as in binary search. 

i j 

Finding the value of k for given i'and j is no longer as 
simple, however; in the actual method used, two arrays of 
indices called Last and Next are precalculated. If U is 

less than P the following comparison is made with P ; 

k Last(k) 

if U is greater than P then P is used for the 

k Next (k) 

succeeding comparison. Termination of the search is 
indicated when Last (k) = 0 (choose event k) or Next(k) = 0 

(choose event k-s- 1 ) . 

The binary search scheme used here is closely 
related to Huffman (minimum length) codes in information 

theory and is essentially similar to Knuth's optimum binary 
search tree algorithm [15]. It can be shown [15] that 
between H(P) and H(P)+2 comparisons are reguired for this 
method, where H(P) is Shannon’s entropy function 



H (P) 



N 

.z,p. 
1=1 1 



log (1/p. ) . 

2 l 



For example, for k=0. 1 it was found that M = 2N + 2 = 164 
and that less than 6.4 comparisons were needed on the 
average to select an event. 



2 • R ec t angle Method 

If the binary search scheme selects rectangle i and 
if V is a second uniform (0,1) random variable independent 
of U, then 
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Z = x + Vh 
i i 

is a random number selected from the i-th rectangle. 



3 • Rejection Method for W ed qes 

Knuth [14] and Marsaglia [23] discuss in some 
detail a method for sampling a random variable whose density 
function is almost linear. This method is used without 
modification to sample from the wedges. 

The values of a and b were found for each wedge 

i i 

so that the gamma density function f ( x) lies within the 
strip indicated in Figure 2 , that is so that 



a 

i 



Since f(x) 



b (x - x )/h < f(x) < b -b(x-x)/h 

i i i i i i i 



for x < x < x 

i i+ 1 



is concave upwards for all x and since the 




Figure 2. Sampling from the i-th wedge. 



algorithm is fastest for a close 

i 

taken to be the tangent to the 



to b , the lower line is 
i 

density function curve 
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parallel to the chord from (x , y ) to (x_ f tJ . The 

i+ 1 i 1 1 

tangent can be found by New t on-Raphson iteration. 

The wedge algorithm is then: 

(1) Generate two uniform random numbers, U and V; 
if G > V then exchange U and V. 

(2) Set Z to x Uh . 

i i. 

(3) If V < r = (a -y ) / (b -y ) go to step (5). 

i i i i i 

(4) If V > U + f (Z)/b go back to step (1) . 

i 

(5) Z is the desired variate; stop. 

When a and b are close together then r will be 
i i i 

nearly one and the algorithm will terminate as a result of 

step (3) most of the time. For roost of the wedges in the 

gamma generator, values of r are greater than .65 sc that 

i 

f (x) needs tc be evaluated (step (4) ) less than 1/3 of the 

time. 

Knuth [14] gives a proof that this method works 

properl y . 

B. APPROXIMATIONS FOR THE INVERSE GAMMA CDF 

Any continuous random variable can be sampled if its 

-1 

inverse distribution function F can be found, for if U is 

-1 

a uniform (0,1) variate then Z = F (U) is a variate from 

the desired distribution. Unfortunately, there is no simple 
inverse for the gamma distribution function with shape 
parameter less than one; Phillips and Beightler [26] 
attempted a rational approximation for the inverse with 
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little success (see Subsection III- E. 2) . In the present 
method, two different inverse approximations are used for 
the head and tail regions- No inverse is needed elsewhere 
since the wedge sampling technique described above uses 
f (o) , the gamma density function, and not the gamma inverse. 



1. Low Values of Z 



a. Method 

The series expression for the incomplete gamma 
function (ganma CDF) can be expanded as follows: 

F (x) = g (k; x) 

co n n+k 

(3-2) F(x) = 1 2 JMI x 

TW n=0 nTln+TT - " 

k 2 

= X [1- x + x 

rW 7. 2+1 'IXZT7T 

If Z = F ( x) then 

k 2 

kF(k)Z = x [1- k x+ k x - ...] 

Ic+T 21kT2T J 

1/k 2 1/k 

[T(k+1)Z] = x [1 - k x + k x - ...] 



When k^_x « 1 the first-order approximation to 
R + i 

the inverse, x , can be used. This corresponds to the case 

1 

when the variable is to be sampled from that part cf the 
head region close to zero. 



x = [ r (k+1) Z ] 

. -1 

= F (Z) 



1/k 
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When the k x term cannot be ignored, 
Tc+T 

approximation is obtained by solving 

2 1/k 

x = x[1-kx+0(x)] 

1 Tc+T 



a better 



= x [ 1 - 1 • k x + 0 ( x ) ] 
Tc T,T T 

2 3 

=x- 1 x +0(x) 



for the second-order approximation, The result is 



x 2 = (1/2) r (k+1) - /TTc+Tf 2 - ] 

-1 

= F (Z) . 



A more convenient form of x for computation is 

2 

x 2 = ( 2x ^ ) / [1 + /T--^ 7 WTy ]. 
b. Accuracy 

The value of P (probability of the head 

0 

region) is determined empirically for the gamma generator by 
the formula 

-5 2 

P = 10 / k 

0 

so that the accuracy of the approximation to the inverse 
incomplete gamma function was investigated for values of Z 

in the range zero to P for various values of k. It was 

0 

found that the approximation improved as Z — > 0 and as 
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-9 

k — > 0; in no case did the relative error exceed 10 in 



-7 

double precision or 10 in single precision. 



2. V al ues of Z Near One 
a. Method 

Newton-Kaphson iteration is used for 
approximating the inverse gamma distribution function for 
values of Z near one; that is the equation 

h (x) = g(k;x) - Z = 0 

is solved for x by using Newton's formula 

x = x - h (x ) /h 1 (x ) . 
n + 1 n n n 

Six-digit convergence can be acheived after four or five 

iterations of this method starting with x = 1.0. 

1 

Application of the method requires evaluation 
of h(x) and h* (x) . The latter is just the gamma density 
function 

k-1 -x 

h * (x) = f (x) = 1 x e , 

while the former can be found by summing the series (3-2). 
Since x is likely to be very large when Z is near one, 
however, as many as 30 terms of the series must be summed 
for convergence with a resultant loss of speed and the 
accumulation of serious round-off errors.' For this reason, 
g(k;x) is evaluated by a continued fraction approximation 
[ 1 ] which has the added advantage of producing the value of 
1 - g(k;x). Since the method is applied only to values of Z 
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in the range (0.999,1-000) the ability to deal with 1-Z 
leads directly to a net gain in significant digits. 

b. Accuracy 

It was difficult to investigate the accuracy of 
this approximation in the same way as for the previous 
method because of the loss in significant digits in 1 - Z as 
Z — > 1. A series of values of y in the range .00001 to 
.00100 was defined; for each value the approximation was 

-1 

used to find x = F (1-y) and then z = F(x) was obtained by 

the continued fraction approximation. The relative error 
between z and 1-y was in most cases zero and in no case 

-7 

greater than h°10 

C. DESCRIPTION OF COMPUTER PROGRAM 
1 * S truct ure 

The gamma generator was implemented as a 

FORTRAN- callable subroutine with each call returning a 
single floating point (type REAL) gamma variate. One 
version was written in the IBM implementation of FORTRAN IV 
while another was written in a combination of FORTRAN and 
IBM 360/Assembler F. Copies of listings for both pregrams 
are on pages 77-04 and 85-88. 

The basic approach in both generators is to divide 
the subroutine into two parts: a relatively lengthy (and 

slow) initialization section which calculates tables and 
constants and a short, fast section which is called 

repeatedly to handle the actual generation process. The 
initialization part of the program was given the name GMINIT 
(for GaMma INITi alization) while the generator was called 
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GAMA. Calling sequences and arguments are discussed in the 
following Subsection. 

Gill NIT accepts values of k in the range 0.05 < k < 
1.00. It produces the cumulative vector PROB (analogous to 
P in Subsection III.A.1); before being placed into PROE the 
component probabilities are ordered with the smallest first 
to minimize round-off errors. GMINIT also finds the arrays 
NEXT and LAST for use in the binary search method as well as 
the arrays X,H,R and B for sampling from the rectangles and 
wedges. Finally, GMINIT calculates constants for the head 
and tail approximations. Sample output from GMINIT is on 
pages 70-72.. 

Two further subroutines were also written: 
IGAM (K,X) which evaluates the incomplete gamma function 
g (k;x) and INVGAM (K,Z) which uses Ne w ton- Raphson iteration 
to solve the equation g(k;x) = 1-Z for x. INVGAM was also 
programmed in both FORTRAN and Assembly language. Listings 
for both IGAM and INVGAM are on pages 89-90 and 91-95, 
respectively. 

In addition, a uniform random number generator is 
required. Since an average of about 2.5 uniform numbers are 
required for each gamma variate, timing characteristics of 
the uniform generator are a significant factor in the 
performance of the gamma generator. The particular 
generator used is called RANDOM and was developed by Lewis, 
Goodman and Miller [19]. RANDOM is very fast; it is capable 
of generating a uniform deviate in about 15 microseconds. 

The generating routine, GAMA, carries out binary 
search with a single uniform random number and then samples 
from the proper sub-distribution using one or more 
additional independent uniform variates. It also performs 
scaling on the generated variate by multiplying it by a 
scale factor BETA which is an input parameter to GMINIT. 

For small values of k (less than 0.1) there is some 
probability of a fixed point underflow error, i.e. that the 
generated variate will be less than the smallest possible 
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-IQ 

floating point number (about 10 for the IBM 360) . The 

FORTRAN standard fixup for this error is to return a result 
of zero , which might be satisfactory for some applications; 
it could create serious problems, however, if, for example, 
a log transformation of the gamma variates is required. If 
a large enough scale factor is used, modification of GAMA to 
perform prescaling before applying the low value 
approximation would reduce the problem. 

2 - Implementati on s 

a. FORTRAN IV Version 

In the FORTRAN version of the generator the 
instruction sequence 

CALL G MI NIT (K, BETA) 

+ m w 

CALL GAMA (Z, IX) 

results in the initialization of the generator for shape 
parameter K and scale parameter \ = 1/BETA and the 
assignment of a gamma random value to Z. IX is an integer 
seed for the random number generator RANDOM and should not 
normally be modified in the main program. The same seguence 
of gamma deviates can be reproduced at any time, however, by 
re-initializing IX to its original value. 

IBM FORTRAN IV has several features that are 
not in American National Standard FORTRAN and which may not 
be implemented in other compilers. Two of these are used in 
the generator and its subroutines: the built-in function 

GAMMA (X) which computes the gamma function of its argument 
and the ENTRY statement which allows a subroutine to have 
multiple entry points. The ENTRY option is used to cut down 
on overhead in differentiating which of the two parts 
(GMINIT or GAMA) of the subroutine is being called; it is 
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equivalent to having two independent subroutines sharing a 
COMMON block- Besides modification of these areas, 
conversion to another computer would require adjustment of 
convergence test values in INVGAM, T.GAH and GMINIT. 

b. Assembler Version 

The Assembler version of the generator was 
written to save execution' time in GAMA; the instruction 
seg uence 

CALL START (K, BETA, IX) 

Z = GAMA (0) 

has the same effect as the previous sequence. Linkage 
overhead in GAMA has been substantially reduced, however, 
since no arguments need be passed (the zero in the 
instruction above is a dummy argument and is required by 
FORTRAN conventions). The capability of re-initializing the 
seed for RANDOM was not retained in this version. 

START is a dummy section of GAMA which calls on 
a modified version of GMINIT. Since execution time was a 
secondary consideration for the initialization routine 
(which was normally called only once) , GMINIT was 
maintained in FORTRAN for ease of programming. The only 
modifications to GMINIT in this case are in the calling 
sequence and in the values produced for the arrays NEXT and 
LAST, which were changed to allow for easier search. 

3 . Tim ing and Core Requirements 



The generator was written with the sole aim of 
achieving the fastest possible execution time without regard 
for storage requirements. The execution times obtained were 
highly dependent on k; for example, more rectangles are 
needed for lower values of k so that the binary search 
method occupies more time in these cases. Samples of 10,000 
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variates were run on the Naval Postgraduate School IBM 
360/67 computer and the generation times recorded for 
several different values of k; the average times per variate 
(in microseconds) were as follows: 



k 


FORTRAN IV Version 




Asse mbler 


Johnk * s 


G 


Compiler 


H Compiler 




Version 


Method 


0.10 


430 


327 




201 


824 


0.25 


385 


310 




189 


864 


0.50 


350 


279 




186 


889 


0.75 


339 


273 




178 


875 


0.90 


325 


270 




179 


80 3 


Av erage 


366 


29 1 




1 87 


851 


By 


way of 


comparison. 


an 


implementation of 


Marsaglia* s 


normal generator on 


the 


same computer 


requires 


about 60 mi 


crosecon ds 


per variate 


while a 


Marsaglia 


exponential 


generator 


takes 70 


microseconds per 


variate. 



Johnk's method is discussed in Subsection III.E.3; it is 
included here for purposes of comparison. 

Although this method of evaluating run times is 
easy to implement and gives a realistic estimate of what can 
be expected in practice, it is somewhat misleading since 
competition for system resources between jobs in a 
multiprogrammed computer causes considerable variability in 
observed run times for a given job. The observations above 
were recorded when the system was relatively idle; on a busy 
afternoon these times may increase by as much as 15 per 
cent. 

The times to complete the set-up phase in GillNIT 
were also highly dependent on k; the following run time 
values vrere observed: 



k 


Set-up Time 
(seconds) 


0. 10 


.893 


0.25 


.348 


0. 50 


.212 


0. 75 


.160 
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Since a much larger number of rectangles is required for 
decomposition when k is small, it can be seen that the 
additional computation time for set-up in these cases can be 
quite costly. In fact, enough time is used to generate 
between 900 and 4500 variates, depending upon the value of 
k. 

Both the Assembler and FORTRAN versions of the 
generator reguire excessive ' core storage, at least in 
comparison with other methods for generating random variates 
(see, for example, Ahrens and Dieter [ 2 ]) . At the Naval 
Postgraduate School installation, where the minimum job core 
allotment is 100,000 bytes (25,000 words) , the size of the 
program was not a restriction, although it could be at 
another installation. Both versions of the subroutine 
require, in addition to space for their own code and the 
tables, core storage for the FORTRAN built-in functions EXP, 
DEXP, ALOG, DLOG , GAMMA, DGAMH A and SQRT, as well as the 
user subroutines IGAM, INVGAM and RANDOM. Core requirements 
(in bytes) for the various subroutines are as follows: 

FORTRAN IV Version Assembler 

G Compiler II Compiler Version 



GAMA 


13206 


12400 


5664 


GMINIT 






6860 


INVGAM 


1070 


866 


684 


IGAM 


740 


574 


740 


RANDOM 


320 


320 


320 


Functions 


5056 


5 056 


5056 


Total 


2039 2 


19216 


19324 



Further improvements in both space requirements and 
run time may still be possible through the use of in-line 
coding instead of subroutine calls. These modifications 
would have the further result of making the generator 
independent of vendor supplied mathematical subroutines. 
Ahrens and Dieter [2] point out some' advantages of this 
approach in their improvement of Marsaglia's normal 
generator. 
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D. STATISTICAL TESTING 



Having thoroughly investigated the numerical accuracy of 
the various approximations used in the generator, it was 
then necessary to test the generated variates for 
"statistical accuracy". Formally, the null hypothesis 

H : The generated deviates z ,z were 

0 12 

sampled from a unit gamma distribution 

with shape parameter k. 

was to be tested against the alternative 

H : The generated deviates were not sampled 

1 

from the unit gamma distribution. 

Of the wide variety of standard goodness-of-f it tests, two 

were chosen for testing H : the chi-square test and the 

' 0 

Anderson-Darling test. In addition, the reproductive 
property of the gamma distribution (that is, that the sum of 

two gamma random variables with shape parameters k and k 

1 2 

and equal scale parameters is also a gamma random variable 

with shape parameter k +k ) was exploited in a test on the 

1 2 

sums of several variates. 

1 • Chi^squar e Test 

A modification of the chi-sguare test suggested by 
Naylor [25] for uniform random number generators was carried 
out. For this test, a set of n quantiles of the object 
distribution (i.e., the gamma distribution) is first found; 
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the i-th quantile r is the solution to the equation 

i. 



F (r ) = i / n . 

i 



Quantiles 
of size N 

v ar iates 

recorded. 



for the gamma distribution are unique. A sample 
variates is then generated and the number of 

f falling into each interval r < z < r is 
i i-1 i 

The statistic 



2 

X 

1 



n 

n Z (f 
U i=1 i 



2 

N ) 
n 



2 

is then found for the sample; has approximately a 

chi-square distribution with n-1 degrees of freedom. This 
procedure is then repeated M times and a second statistic 
2 

X^ is calculated in an analogous manner, using m quantiles 

of the chi-sguare distribution with n~ 1 degrees of freedom. 

2 

The value of X is then used to accept or reiect H . 

2 0 

Values chosen for the actual test were m = n = 10, 

N = 1000 and M = 100. Thus, each test involved generation 

of 100,000 variates. Tabulated results for various values 

of k are as follows: 





2 




2 




k 


X 


Chi-square 


Minimum X 


Maximum 




2 


Level 


1 




06 


4.8 


. 149 


2.24 


27.28 


10 


4.6 


. 132 


1.80 


22.20 


15 


11.2 


.738 


1.90 


23.74 


50 


4.4 


.117 


1 .66 


21.58 


90 


11.4 


.751 


1.66 


25.68 


95 


8.8 


.544 


1.90 


21.98 
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The column headed "Chi-square Level" gives the quantile of 
the chi-square distribution with nine degrees of freedom at 

2 

the observed value of X . Eased on these results, H could 

2 0 

not be rejected. 

There are some serious drawbacks to the chi-square 
test, not the least of which is the fact that the statistic 
depends strongly on the specific quantiles taken and is not 
an invariant for a particular sample. Nevertheless, the 
test proved to be quite sensitive even to minor variations 
in the generator. For example, it readily detected the 
statistical effect of the omission of a single wedge (whose 
probability was .002) by the binary search procedure, 
producing composite scores greater than 50 in this case. As 
is indicated in the sequel (Subsection III.E.2) , the test 
was also used to reject the Phillips generator, which 
acheived scores greater than 800 for small values of k. 

2. Anderso n- Darling Test 



The test of good ness- of-f i t proposed by Anderson 
and Darling [3] is a distribution-free test in that instead 



of the observed values of 
val ues 



the variate x ,x , 
1 2 



x the 
n 



u = F (x ) 
i i 



are used. If F(x) is in fact the distribution function of 



the x • s (that is, if the null hypothesis is true) , then the 
i 



u 1 s are a sample of size n from 
i 



distribution. If the u *s are now ordered 

i 



a 



so 



uniform (0,1) 



that u 

( 1 ) 



< 
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u < ... < u , then the Anderson-Darling statistic 
(2) (n) 

2 n 



W = - n - 

n 


1 .2 r 

n i = 1 


(2i - 1) In 


u 

(i) 








+ {2 (n- 


i) +1} In 


(1 - 


u . ) 
(i) 


] 


can be computed. The 


null hypothesis can 


then 


be 


accepted 


2 

if W is not too 


large. 


Lewis [17] 


has 


determined the 



n 

2 

distribution of W for various values of n. 

n 

It should be mentioned that the Anderson-Darling 
statistic gives the most weight to observations from the 

tails of the distribution, where u is very close to zero 

(i) 

or one. Thus, the test was ideal for the gamma generator, 
since the rectangle-wedge technique had been successfully 
applied in a wide number of other cases and the main 
uncertainty for the generator lay in the statistical 
properties of the head and tail approximations. 

Accordingly, Anderson-Darling statistics were found 
for 100 samples of 100 variates each. The results of Lewis 
[17] were then used to find ten guantiJ.es of the 

2 

distribution of W so that a chi-square statistic could 

100 

be computed for the sample of 100 Anderson-Darling 



statistics. The results were as 


follows: 




k 


Chi-square 

Statistic 


C hi- square 
Level 


2 

Minimum W 

n 


Maximum W 


. 10 


4.269 


. 107 


. 178 


3.500 


. 50 


4.895 


. 157 


. 175 


5.071 


. 90 


4.306 


.110 


.223 


5.582 



4 5 



Once again 



test. 



it was decided to accept II on the basis of the 

0 



3 • " Reprod u cti ve 11 Pro joe r t_y Test 

It is well-known (Cramer [8]) that the sum of two 
independent gamma random variables with identical scale 

parameters and respective shape parameters k and k also 

12 

has a gamma distribution with shape parameter k +k . Since 

12 

a gamma random variable with k=1.0 is an exponential 
variable, a possible test for is to investigate whether 
the statistic 



V 

n 



n 

z 

i= 1 l 



r 



has a unit exponential distribution, where the z 's are 

i 

generated by GAMA with k = 1/n. 

The test was carried out for samples of 100 V { s 

n 

with n = 2,4,6,8,10. For each of four such samples, an 
Anderson-Darling statistic was computed and the maximum 
value of the statistic recorded. The approximate 
significance level of the statistic was determined by linear 
interpolation in the table of the asymptotic distribution of 
2 

W given by lewis [17]. The results are as follows: 
n 
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2 







n 


k 


Maximum W 


A 


ppro ximate 
















100 Signi 


ficance Leve 


1 








2 


.500 


1.577 




. 84 1 










4 


.250 


2.045 




.913 










6 


. 167 


1.031 




.658 










8 


. 125 


1.469 




.816 










10 


. 100 


1.525 




. . 829 










As 


a further test 


of whether 


the statisti 


cs V 

n 


have 


an 


expon 


ential d 


istr ibution 


, the stati 


stic 












V' 


= Minimum (V , V . . 


1 • f V ) 












n 




n ; 1 n ; 2 


n; 1000 






was 


inve 


sti 


gated 


; this was 


done speci 


fically to 


test 


the 


statistical 


effect of the low value approximation 


-1 

to F 


(Z) . 


It 


is we 


li- 


known 


that the 


minimum of 


m unit e 


xponcn 


tial 


ran 


dom v 


ar i 


ables 


is an expo 


nential random variable 


with 


mean 


1/m 


. Th 


us, 


samples of 100 


observation 


s of V s were 
n 


generated 


for 


n = 2 , 


4 ,8 


and 


an Anderson 


-Darling st 


atistic comp 


u ted u 


sing 


the 


expo 


nen 


tial 


distribu tio 


n function 


with mean 


001. 


The 



results are as follows: 



2 



n 


k 


W 

100 


Ar 

Signia 


>proximate 
ricance Level 




2 


0.500 


.903 




.588 




4 


0. 250 


. 431 




.182 


t 


8 


0. 125 


2.287 




.936 





As can be seen, H is readily accepted based on these tests. 






COMPARISON WITH OTHER METHODS 



E. 



1 • Na y l or * s M et hod 

The method of Naylor [25] for producing gamma 
variates with non- in tegral shape parameter was not seriously 
considered for the simulation since it is not defined for 
k < 1,0. The method involves sampling a mixture of two 

gamma random variables with integral shape parameters L kJ 
and l]cj + 1 (the notation L k J denotes the truncated integer 

part of k) . Berman has shown [28] that the method is 

satisfactory for small simulations when k>5 but his results 
demonstrate that the method is not statistically acceptable 
for sample sizes on the order of 100,000 for any value of k. 

The timing results for the Naylor method obtained 
by Berman [28] for an IBM 360/65 computer suggest that a 

better way to generate non-integral gamma random numbers is 
to add an integral gamma variate with shape parameter L k J 
and a variate from GAMA with shape parameter k - L k-* . For 
k>10 this modification adds less than 20 per cent to the 
execution time but renders the result statistically exact. 

2 . P hi l lips* s Met hod 

Phillips and Beightler [26] proposed the use of a 
rational approximation to the inverse gamma distribution 
function for generating gamma variates; three different 
approximations were developed for 0 < k < 2.0, 2.0 < k < 5.0 
and k > 5.0. The technique was extensively tested by 

Phillips using the Kolmogorov-Smirno v goodness of fit test; 
unfortunately, this test is notoriously insensitive to 
deviations in the tails of the distribution. The chi-square 
goodness of fit test was therefore applied to the method for 
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k < 1.0 and it was found that the hypothesis that the 
Phillips numbers were sampled from the gamma distribution 
could be rejected at any desired confidence level 
(chi-square scores in excess of 800 on 9 degrees of freedom 
were observed). The method was not investigated for k>1.0. 
Furthermore, this method required substantially more 
execution time than the method developed here (on the order 
of 600 microseconds per variate). 

3 . 02li£]S_l2 Me th od 

Johnk's method [13] was the only known 
theoretically exact means of generating gamma variates with 
non-integral shape parameter. It is a rejection technique, 
like the wedge method of of Subsection III. A. 3. For 
generating gamma variates with shape parameter k less than 
1.0 the method can be set forth as follows: 

(1) Generate two independent uniform (0,1) random 
variates U and V. 

1/k 1/ ( 1- k) 

(2) Set X = U and Y = V 

(3) If X + Y > 1.0 go back to (1). 

(4) Generate a new uniform variate R and set 
W = In R. 

(5) Z = W * X / (X+Y) is a unit gamma variate. Scale 
Z if required and stop. 

Translations of Johnk's original proofs of the correctness 
of this methcd are contained in the Phillips and Beightler 
[26] and Fox [9] papers. Essentially the method generates a 
beta variate and then uses the Lukacs-Laha results [16,21] 
to produce a gamma variate. Johnk shows that when k is less 
than 1.0 the rejection in step (3) will occur between 1.00 
and 1.27 times per variate on the average. 

Johnk's method is not satisfactory for large-scale 
simulation use because it requires excessive time; see 



Subsection III.C. 3 for the observed times for a FORTRAN 
implementation of the method. The time can be improved by 
recognizing that if U is a uniform (0,1) variate then -In U 
is a unit exponential variate. The time to generate U and 
evaluate its natural logarithm is about 130 microseconds on 
the IBM 360/67 computer while an assembler implementation of 
Marsaglia's exponential algorithm can deliver an exponential 
variate in about 70 microseconds. The improved algorithm is 
the n 



(1) Generate two independent exponential variates, R 

and S . 



-R/k -S/(1-k) 

(2) Set X = e and Y = e 

(3) If X + Y > 1 go back to (1). 

(4) Generate a new exponential variate T. 

(5) Deliver Z = T • X /(X+Y). Stop. 



Even with these improvements, however, the method still 
requires ever 500 microseconds per variate in an assembler 
version. Thus, allowing for GAHA's substantial set-up time, 
it is faster to use the Johnk method for less than 250 
variates (k=0.75) to 1,350 variates (k = 0.1), depending on 
k, but the method is too slow for large-scale simulation 



(100,000 


or 


more 


variat es) . 


The 


method remains very 


at tractive 


f cr 


use 


when time is 


not so 


critical, however. 


since it 


is 


ver y 


easy to 


program 


and requires minimal 


storage. 
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IV. COMPUTER SIMULATIO N EXPERIMEN T 



A. DESCRIPTION OF METHOD 

The basic purpose of the simulation experiment was to 
investigate the distribution of Y for small values of J. 

ij 

Different si mutations were run for shape parameters k of .1, 
.25, .5 and ..75 corresponding to coefficients of variation 

2 

C (X) of 10 f 4, 2 and 1.33. In addition, a simulation using 

an implementation of the Marsaglia exponential generator was 
run for the case k=1.0. 

1 . Simulati on 

Ten thousand samples of J gamma variates each were 

generated for J - 10,30,50,100 and the Y statistic computed 

J 

for every sample. The first four sample moments were found 

and then the entire sample of 10,000 Y 1 s was ordered in 

J 

order to estimate the quantiles of the distribution (see 
Subsection IV. B. 2). This procedure produced 27 statistics: 
four sample moment estimates and 23 quantile estimates. The 
entire procedure was then repeated ten times and the mean 
and variance of each of the 27 estimates found over the ten 
samples of each. 

Obviously a great many gamma variates were reguired 
to carry out the simulation, nearly 20 million for each 

value of k. The size of the simulation graphically 

demonstrates the need for speed in the gamma generator. 
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Even with the relatively fast performance of the 
Assembly- coded generator, the total IBM 360/67 central 
processor time needed to run the simulation was about 60 
minutes for each k value. 

2. FORTRAN IV P rog ram 

The simulation itself was written in FORTRAN IV; a 
program listing is included on pages 96-98. The program and 
subroutines required nearly 89,000 bytes of core at run 
time, 40,000 bytes of which were needed for the 
10,000-element array used for quantile estimation. A sample 
of output frcm the program for J=30 is included on pages 
73-76. 

B. ESTIMATION OF PARAMETERS 
1 • bistri bu t ion Mo men ts 



If Y represents the i-th observation of Y in 
J; i J 



the simulation, then the r-th moment m of Y can be 

r J 

estimated by 



m 

r 



N 

1 1 (Y ) 

H i=1 J;i 



An unbiased estimator of the variance is then 



_2 

s 



N 

H=‘ 



(m - m ) 



In the simulation, estimates of the first four moments 
(r=1,2,3,4) were obtained. In addition, the sample variance 
and the coefficients of skewness and kurtosis were 
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calculated from the moment estimates. 



2. Quantiles 

a. Quantile Estimation Considerations 

A quantile x of a distribution F(x) is defined 
a 

as the solution to the equation 

F (x ) = a (0 < a < 1 ) . 

a 

Quantiles are unique for most continuous distributions but 
are not, in general, for discrete distributions. In the 

simulation, quantiles of Y were estimated for a = .001, 

J 

.002, . 005, .010, .020, .025, . 050 , . 1 ( . 1 ) . 9 , . 950 , .975, 

.980, .990, .995, .998 and .999. 

There are two principal methods for quantile 
estimation in simulations: the order statistic estimator 

and the Robbins-Honro stochastic approximation technique. 
(For a more complete discussion, see Goodman, Lewis and 
Robbins [11 ]. ) The order statistic method was chosen for 
this simulation because of the ease of implementation. Let 

the H observations Y , Y , ... , Y of Y be ordered 

J ; 1 J ; 2 J ; N J 



so that Y 


>H 

VI 


< ... 


< Y 




; then Y 


is 


(J;i) 


(J;2) 




(J; n) 




(J ;i) 


called the i-th 


order 


sta tistic 


of 


the 


sampl e 


and an 



estimate of the a-guantile y of the distribution of Y is 

a J 



y = Y 

a (J ; L aU J ) 
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b. Sorting Random Data 



In order to estimate th 



the distribution of Y it is necessary 

J 



samples of Y . 

J 

Subsection) also 
previously, a sain 
est i ui at icn. 

Sorti 

considerable overh 
t he sorting met 

bewildering variet 
purposes (Knuth [1 
The particular tec 
to Singleton [27]; 
capable of orderin 
numbers in less 
therefore, less th 
spent in sorting t 



Considerations of 



call for large sa 
pie size of 10,000 w 



ng such a large 
ead into the simulat 
hod is carefully 
y of known sorting 
5] devotes an entire 
hnique selected for 
a FORTRAN implement 
g a sample of 10,000 
than four seconds 
an five per cent of 
he samples. 



3. Errors and Bias 



e extreme quantiles of 


to deal with large 


bia s 


(see the next 


inples. 


As mentioned 


as c ho 


sen for quantile 


sample 


can introduce 


ion. 


however, unless 



chosen. There is a 
methods for various 
text to the subject) . 
the simulation is due 
ation of the method is 
normalized type REAL 
. In the simulation, 
the execution time was 



The parameter estimates may not correspond to the 
true parameter values for three reasons: statistical 
variability, computation errors and estimator bias. 
Statistical variability is an inescapable reflection of the 
fact that a simulation is not a deterministic process, as is 
the summing cf an infinite series. Thus, the probability 
that the precise expected value of an estimator will be 
observed is very low and may be zero for a continuous 
distribution. The law of large numbers, however, guarantees 
that the estimator will converge closer to its theoretical 
value as more replications of the simulation are performed. 

Computation error is a reflection of digital 
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computer round-off errors; that is, for every computer 
system there is a floating-point constant e such that when 

I d | < | e | 

1.0 + d = 1.0. 

-7 

For the IBM 360 system, e = 9.53 © 10 . Thus, if very many 

data elements are to be summed (as in determining the moment 
estimators) there may be some loss of significance as the 
accumulated sum increases. The problem will be even more 
acute if very skewed data are to be added, as in the 

determination of Y and Y 

1J 2 J 

Bound-off errors in Y and Y , however, tend to 

1 J 2 J 

cancel out when Y is calculated. On the other hand, in the 

J 

simulation the observed moment estimators 1 were 

r 

consistently below their theoretical values, thus reflecting 

considerable round-off error. In fact, m was less than m 

r r 

in 17 observations out of 20 for r=1 and in 14 of 20 for 
r=2 . Computing the moment estimators using double precision 

-16 

arithmetic (thus decreasing the value of e to 2.22 ® 10 

and therefore greatly reducing round-off errors) resulted in 

values of m more in agreement with theory, 
r 

Estimator bias is a theoretical property of the 
specific estimators used in a simulation. In this case, the 

ni are unbiased estimators of the m but 
r r 



E[y 3 
a 



y + 

a 



0 ( 1 ) , 
u 
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where N is the sample size (Goodman, Lewis and Bobbins 

[11])- Thus y is a biased estimator. No attempt was made 
a 

to evaluate the extent of the bias of y , but it is clear 

a 

that the simulation results may not be accurate for more 
than two or three decimal places. 

C- SIMULATION RESULTS 



1 . Agree me nt wi th Theory 



For each simulation run, a t-statistic was computed 
using the sample mean and variance. Since the sample mean 
is the sum of 100,000 identically distributed random 

variables (the observations of Y ) , it is safe to assume 

that its distribution is very nearly normal and that 



t = <m - EC Yj 3) / (s / /n) 



where 



_2 

s 



n 

n^ 




m 



2 



1 



) 



has a t distribution with n degrees of freedom. When 
large (100,000 in this case) then the t distribution 
approximated by the standard normal distribution, 
statistics for the various runs are summarized 
following table: 



n is 
can be 
The t 
in the 
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J 

k 


10 


30 


50 


100 


0. 10 


-0. 0688 


-2.2718 


-0.564 1 


-0.9348 


0. 25 


1 . 0851 


-2.74 54 


-2.4309 


-1.9453 


0.50 


0- 1989 


-3. 6603 


-2.5376 


1.2218 


0.75 


-0.1440 


-1 .8950 


-2.6556 


-1.2765 


• 

O 

O 


-0.5925 


-3.4506 


-1.2623 


NA 



Since the .01 quantile of the standard normal distribution 
is -2-326, it can be readily seen that the above results do 
not agree with theory. 

As mentioned above, the reason for this discrepancy 
lies in accumulating round-off errors in calculating the 
moment estimators. The simulations for J=30 were thus rerun 
with the moment estimators computed in double precision. No 
change in the quantile estimates was observed, but the new t 
statistics were 

k 0.10 0.25 0.50 0.75 

t statistic -0.9206 -0.8101 -1.0592 1.3176 

Thus, it is concluded that the simulation results agree with 
theory as far as the observed moments are concerned. 

2 . Tabulat io n of the Distr i but ion of Y 

J 



In addition to the actual values for Y , the 

J 

simulation program produced normalized values for the 
quantile estimates 



n = (y - u) / a , 
a a 



w here 
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u = Ef Y ] 
J 



J+ 1 



and 



c 



z 



Var£ Y ] 
J 




j+1 - 
TcJ+T 



The normalized statistics are tabulated here in 

show the asymptotic convergence of the distribution 

a normal distribution; quantiles of the standa 
distribution are also listed* 



order to 

of Y to 
J 

d normal 
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Normal 

Quantile 



Alpha J =JK) J=30 Jf 50 J = 1 00 



.001 


-2. 2027 


-2. 7519 


-2. 8819 


-2. 9941 


-3.090 


.002 


-2. 1920 


-2.6114 


-2.7360 


-2.7978 


-2.878 


.005 


-2. 1473 


-2. 4244 


-2.4940 


-2.5415 


-2.576 


.010 


-2.0788 


-2. 2439 


-2.2911 


-2.3166 


-2. 326 


.020 


-1.9445 


-2. 0244 


-2.0507 


-2. 0549 


-2.054 


. 025 


-1.8738 


-1.9473 


-1 .9655 


-1.9672 


-1.960 


.0 50 


-1.6546 


-1. 6654 


-1.6675 


-1.6547 


-1.645 


. 100 


- 1. 3444 


-1.3222 


-1 .3107 


-1.2983 


-1. 282 


.200 


-0.9230 


-0. 8844 


-0. 8697 


-0. 8563 


-0.842 


. 300 


-0.5885 


-0. 5541 


-0.5441 


-0.5364 


-0. 524 


.4 00 


-0. 2773 


-0. 2746 


-0.2608 


-0. 2613 


-0.253 


. 500 


-0.0012 


-0. 0051 


0.0023 


-0.0029 


0 


.600 


0.2802 


0. 2687 


0.2685 


0.2604 


0.253 


.700 


0.5923 


0.5535 


0.5535 


0.5382 


0. 524 


.800 


0.9257 


0.8783 


0.8732 


0. 8566 


0.842 


.900 


1.3426 


1.3177 


1 .3042 


1 . 2958 


1. 282 


.950 


1.6541 


1. 6553 


1.6621 


1.6577 


1.645 


.975 


1. 8754 


1.9228 


1 .9537 


1 . 9627 


1.960 


.980 


1 . 9437 


2.0036 


2. 0467 


2. 0523 


2.054 


.9 90 


2.0777 


2.2178 


2.2859 


2.3096 


2. 326 


.9 95 


2. 1498 


2. 4068 


2. 5013 


2.5484 


2.576 


.998 


2.1901 


2. 6027 


2.7520 


2. 8258 


2.878 


.9 99 


2.2009 


2. 7 274 


2.9476 


3.0085 


3.090 


u 


5.5 


15.5 


25.5 


50. 5 




a 


2. 0310 


4. 3277 


5.8914 


8.7034 




t-'Test 


-0.0688 


-2.2718 


-0.5641 


-0.9348 





for Mean 

Table I. Simulation results for k=0.10. Normalized 

quantiles of the distribution of Y under the null 

J 

hypothesis b=0 (no trend) are tabulated. 
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Al ph a 

.001 
.002 
.005 
.010 
. 020 
.025 
.050 
.100 
. 200 
.300 
.400 
.500 
.600 
.700 
. 800 
. 900 
.950 
.975 
.980 
.9 90 
. 995 
.9 98 
.999 

u 

G 

t-Test 
for Kean 



J=J0 

-2. 6635 
- 2. 5614 
-2. 3792 
-2.2094 
-2.0089 
-1 . 9340 
-1.6581 
-1.3211 
-0. 8811 
-0.5511 
-0.2643 

0.0089 
0.2673 
0.5619 
0.8913 
1.3230 
1.6566 
1.9359 
2.0101 
2.2225 
2. 3901 
2.5484 
2.6496 

5.5 

1.5353 

1.0851 



J=30 

-3.0231 
-2.7949 
-2.5528 
-2. 3081 
-2. 0596 
-1. 9722 
-1.6723 
-1. 3174 
-0. 8663 
-0. 5414 
-0. 2608 
-0.0016 
0.2573 
0. 5365 

0. 8639 

1. 3024 
1.6636 

1. 9642 
2. 0579 

2. 2991 
2.5438 
2.7867 
2.9408 

15.5 

2.9688 

-2.7454 



J=5 0 

-2.9940 
-2.8016 
-2.551 9 
-2.3277 
-2.0612 
-1.9727 
-1 .6639 
-1.2980 
-0. 8604 
-0.5375 
-0.2632 
-0.0070 
0.2568 
0.5311 
0.8589 
1.3015 
1 .6593 
1.9692 
2.0545 
2.3178 
2.5635 
2. 8170 
3.0409 

25. 5 
3.9276 

-2.4309 



J=100 

-3.0984 
-2.8912 
-2.5871 
-2.3405 
-2.0672 
-1.9775 
-1.6715 
-1.2982 
-0.8583 
-0. 5339 
-0. 2585 
0. 0015 
0.2573 
0.5351 
0.8529 
1.2946 
1.6576 
1.9776 
2.0682 
2.3388 
2.5759 
2. 8870 
3.0975 

50.5 

5.6611 

-1.9453 



Normal 

Quantile 

-3. 090 
-2.878 
-2.576 
-2.326 
-2.054 
-1.960 
-1.645 
- 1.282 
-0. 842 
-0.524 
-0. 253 
0 . 

0. 253 
0.524 
0. 842 
1.282 
1. 645 
1.960 
2.054 
2.326 
2.576 
2.878 
3.090 



Table II. Simulation results for k=0.25. Normalized 

quantiles of the distribution of Y under the null 

J 

hypothesis b=0 (no trend) are tabulated. 
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Aljgha 


J = 1 0 


Jf 30 


J = 50 


J^JOO 


Normal 

Quantile 


.001 


-2.8514 


-3. 0592 


-3. 0958 


-3.0890 


-3.090 


.002 


-2. 7088 


-2. 8700 


-2.9104 


-2.9028 


-2.878 


.0 05 


-2.4928 


-2. 5669 


-2.6023 


-2.5900 


-2.576 


.0 10 


-2. 2861 


-2.3410 


-2.3446 


-2.3538 


-2. 326 


.020 


-2.0416 


-2. 0851 


-2.0823 


-2.0925 


-2.054 


.0 25 


-1.9578 


-1.9918 


-1 .9941 


-1.9955 


-1.960 


.050 


-1.6617 


-1. 6761 


-1.6750 


-1.6605 


-1.645 


. 100 


-1. 3028 


-1.3107 


-1.3119 


-1.2956 


- 1. 282 


.200 


-0.8685 


-0. 8723 


-0. 8620 


-0.8468 


-0.842 


. 300 


-0.5420 


-0.5460 


-0.5375 


-0.5201 


-0. 524 


.400 


-0.2661 


-0. 2660 


-0. 2587 


-0. 2471 


-0.253 


.500 


0.0030 


-0.0051 


0.0022 


0.0092 


0 . 


.600 


0.2665 


0. 2592 


0.2596 


0. 2708 


0.253 


.700 


0.5483 


0.5375 


0.5380 


0.5430 


0. 524 


.800 


0. 8727 


0. 8610 


0.8622 


0.8679 


0. 842 


. 900 


1 . 3098 


1.3037 


1.3123 


i . j 1 7z 


't o n n 

i d o zL 


.9 50 


1.6615 


1. 6712 


1.6747 


1.6861 


1.645 


.975 


1. 9576 


1.9847 


1.9851 


2.0026 


1.960 


.980 


2.0483 


2. 0753 


2.0827 


2. 0992 


2.054 


.990 


2.2750 


2.3388 


2.3416 


2.3790 


2.326 


.9 95 


2.4723 


2. 5856 


2.5955 


2.6282 


2.576 


.998 


2.6926 


2. 8816 


2.8901 


2. 9099 


2.878 


.9 99 


2.8498 


3. 0650 


3.0638 


3. 1490 


3.090 


u 


5.5 


15.5 


25.5 


50.5 




a 


1. 1726 


2. 1639 


2.8301 


4.0421 




t-Test 


0. 1989 


-3.6603 


-2.5376 


1.2218 





for Mean 



Table III. Simulation results for k=0.50. Normalized 

quantiles of the distribution of Y under the null 

J 

hypothesis b=0 (no trend) are tabulated. 
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Al_gha 

.001 
.002 
. 005 
.010 
. 020 
.025 
.050 
. 1 00 
.200 
.300 
.400 
.500 
.600 
.700 
. 800 
.9 00 
.950 
.675 
. 980 
.9 90 
.995 
.9 98 
.9 99 

u 

a 

t-Test 
for Mean 



J = 10 

-2. 9330 
-2.7684 
-2.5103 
-2. 2857 
-2. 0410 
-1 . 9544 
-1.6617 
-1 . 3090 
-0. 8670 
-0.5429 
-0.2628 
0.0024 
0.2618 
0.5410 
0. 8707 
1.3153 
1.6657 
1.9598 
2. 0509 
2.2908 
2.5102 
2.7933 
2. 9584 

5.5 

0. 9852 
-0. 1440 



Jf30 

-3. 1428 
-2.9370 
-2.6022 
-2. 3462 
-2. 0911 
-1. 9931 
- 1.6686 
-1.3054 
-0. 8556 
-0. 5275 
-0.2496 
0. 0079 
0.2658 
0. 5380 
0. 8640 
1.3052 
1.6771 
1. 9889 
2.0705 
2.3475 
2.6129 
2. 8948 
3. 1 079 

15.5 

1.7855 

-1. 8950 



J=50 

-3.1428 
-2. 9344 
-2.5863 
-2.3430 
-2.0684 
-1.9815 
-1.6722 
-1.3063 
-0.8629 
-0.5351 
-0.2578 
0.0003 
0.2603 
0.5358 
0.8574 
1.3056 
1.6803 
1.9920 
2.0849 
2. 3639 
2.6022 
2.9219 
3. 0848 

25. 5 
2.3257 

-2.6556 



J=100 

-3.1331 
-2. 9034 
-2.6140 
-2. 3620 
-2.0876 
-1.9907 
-1.6678 
-1.3039 
-0. 8497 
-0.5279 
-0.2495 
0. 0064 
0. 2631 
0.5390 
0. 85 84 
1.3000 
1.6708 
1.9943 
2.0894 
2. 3518 
2.5834 
2.8904 
3. 1 105 

50.5 

3.3112 

-1.2765 



Normal 

Quantile 

-3.090 
-2.878 
-2.576 
-2.326 
-2.054 
-1.960 
-1.645 
-1.282 
-0. 842 
-0. 524 
-0. 253 
0 

0. 253 
0.524 
0. 842 
1.282 
1. 645 
1.960 
2.054 
2.326 

2. 576 
2.878 

3. 090 



Table IV. Simulation results for k=0.75. Normalized 

quantiles of the distribution of Y under the null 

J 

hypothesis b=0 (no trend) are tabulated. 
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J = 100 



Normal 

Quantile 



Al£ha J =10 J=30 J=50 



.001 


-2. 9858 


-3. 0964 


-3.0550 NA 


-3. 090 


.002 


-2.7774 


-2.8721 


-2.8522 


-2. 878 


.0 05 


-2. 5046 


-2.5750 


-2. 5525 


- 2.576 


.010 


-2.2994 


-2.3304 


-2.3006 


-2. 326 


. 020 


-2.0442 


-2. 0496 


-2.0350 


-2.054 


. 025 


-1.9502 


-1 . 9566 


-1.9455 


-1.960 


.050 


-1.6517 


-1. 6464 


-1.6489 


- 1. 645 


. 100 


- 1. 2917 


- 1.2846 


-1.2775 


-1. 282 


.200 


-0. 8543 


-0. 8440 


-0 . 8348 


-0. 842 


. 300 


-0. 5284 


-0.5244 


-0.5172 


-0. 524 


.400 


-0.2540 


-0. 2533 


-0. 2409 


-0.253 


.500 


-0.0009 


-0.0004 


0.0075 


0 


.600 


0.2562 


0. 2566 


0.2607 


0.253 


.700 


0.5286 


0.5238 


0.5317 


0. 524 


.800 


0. 8519 


0. 8422 


0.8443 


0.842 


. 900 


1. 2905 


1.2875 


1 .284 3 


1. 282 


.950 


1.6536 


1.6417 


1 . 6508 


1.645 


.975 


1.9521 


1.9635 


1 .9690 


1.960 


.980 


2.0417 


2. 0588 


2.0541 


2.054 


.9 90 


2. 2855 


2.3152 


2.3165 


2. 326 


.9 95 


2.4933 


2. 5520 


2. 5828 


2. 576 


.998 


2.7571 


2. 8427 


2.8756 


2. 878 


.9 99 


2.9209 


3. 0438 


3.0758 


3.090 


u 


5.5 


15.5 


25.5 




a 


0.8660 


1. 5546 


2.0207 




t-Test 


-0. 5925 


-3.4506 


-1.2623 





for Mean 

Table V. Simulation results for k=1.00. Normalized 

guantiles of the distribution of Y under the null 

J 

hypothesis b=0 (no trend) are tabulated. 
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D. LARGE-SCALE SIMULATION TOOLS 



In programming the simulation it became clear that a set 
of standard data-handl ing subroutines would have been most 
helpful for parameter estimation. The facilities of 
standard simulation languages (e.g., GPSS or CSMP) and data 
reduction programs (SNAP/IEDA or SPSS) are not suited to the 
amount or type of information generated in large-scale 
statistical simulations such as this one. Extensive 
research in the literature was thus required to find 
high-speed techniques for sorting, quantile estimation and 
statistical testing and considerable effort was expended in 
programming the techniques. 

The sophistication of state-of-the-art methods in random 
number generation, for example, can result in substantial 
savings of computer time and core requirements but only at 
the expense of program complexity. Since the specific 
details of the method involved are usually of peripheral- 
importance to the simulation itself, removal of the 
programming burden will result in more resources becoming 
available for the main purpose of the simulation. 

Following Lewis* discussion [24] of a large-scale 

computer-aided statistical package, the following routines 

would have been of value in the Y simulation: 

J 

(1) A fullv-docuraented quantile estimation 
subroutine using the maximum transformation Robbins-Monro 
stochastic approximation technique. 

(2) One or more efficient goodness-of -f i t test 

subroutines: t-test, chi-square test, Anderson-Darling 

test, Kolmogcrov-Smirnov test. 

(3) A sorting routine equal in speed to ACM 
Algorithm 347 [27 ]. 
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moment 



(*0 



est imators 



An efficient subroutine for accumulating 
with minimal round-off error. 
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V. CONCLUSION 



A. APPLICABILITY OF RESULTS 



1 . Gamma Generator 

In view of the wide variety of physical events 
which can be successfully modeled empirically as gamma 
renewal processes, the availability of an exact, high-speed 
gamma variate generator should enable more extensive 
simulations to be conducted in many applications areas. If 
the shape parameter remains constant throughout a given 
application, pre-calculation of GAMA's tables and constants 
will save set-up time in succeeding runs; if not, then the 
trade-off between GAMA and the Johnk technique (Subsection 
III.E.3) should be considered in selecting a generator. 
Nevertheless, it is clear that whenever a large number of 
gamma variates with shape parameter less than one is 
required, GAMA is superior to any other known method for 
generating them. 

2 * T es ts f or Tren d 



variance 
needs t 
squared 
now be 
process 
not . F 
computer 



The results of Section II. D indicate that the 
of the standard B test for the Poisson process 
o be "inflated" by the coefficient of variation 
for use in a gamma renewal process. Thus, it may 
possible to accept the null hypothesis that a given 
is without trend when prior to this result it was 
or example, Lewis and Shedler [20] in analyzing 
page exception data obtained B values which were 
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far too high to accept the null hypothesis (normalized 
values of -2.83, -8.67 and -18.11 were observed; the .01 
quantile is -2.33). However, allowing for the extremely 
high coefficients of variation involved (3.34, 3.27 and 
3.70), the null hypothesis can be accepted in two of the 
three cases. 

Since in practice the actual value of k is seldom 
known, the results obtained in Chapter II suggest that the 
estimated sample coefficient of variation could be used as a 
multiplier of the Poisson variance in applying the E test. 
The simulation results indicate that this procedure may lead 
to slightly higher acceptance levels in the gamma case when 
k is small, but the effect diminishes with increasing J. 

Without investigation of the relative power of the 



Y test against other standard trend tests it may not be 
J 

wise to apply it to an arbitrary process. It would appear 
most useful, in any case, when small values of k are 
involved, for this is precisely the situation which the B 
test does not handle well. 

As for the tabular results. Tables I - V are 
sufficiently accurate for routine application of the test 
based on B, especially since the level of the test is 
arbitrary. Possible inaccuracies in the tabular values 



, if 


the 


valu 


es were used 


to 


of 


the 


Y 


test against 


the 






J 




not 


been 


done. 







The distributions of Y are in all cases observed 

J 

to be approximately symmetric and underdispersed with 
respect to the normal distribution. Convergence to the to 
the normal distribution is rapid, the rate being slowest, as 
expected, for k=0.1. Even there, the normal approximation 
is adequate for J=50 if a level of no higher than 0.990 and 
no lower than 0.010 is required. The approximate rates of 
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to use 



convergence of Y are such that it appears possible 

the normal approximation whenever J is greater than the 
following values: 

k 0. 10 0.25 0.50 0. 75 

J 100 50 30 30 

B. AREAS FOB FURTHER STUDY 

1 . Gamma Gener a to r 

The problem of a general purpose gamma generator (k 
> 1 as well as k < 1) should be investigated in order to 
extend the present work to wider applications areas. A 
straightforward way to do this is to follow Berman's method 
[28] of adding L k-* exponential variates to a gamma variate 
with shape parameter k - L k-» . This approach has two 
distinct draubacks in the present case: 

(1) When k - L k- t is small, GAMA requires much more 
setup and generation time than when 0.5 < k - L k-i < 1.0; 
when k - «-kJ < 0,05, GAMA cannot be used at all. This 
difficulty could, be overcome by using the fact that the 
square of a standard normal random variable has the gamma 
distribution with k=0.5, 1=0.5; thus, if Z is a standard 

2 

normal variate then Y = 2Z is a unit gamma variate with 

shape parameter 0.5. If a Marsaglia normal generator is 
available, Y can be readily generated and added to a series 
of exponential variates to allow GAMA to be called with 
shape parameter in the optimum range. A thorough 
investigation of the time requirements for the various 
options is needed to find the fastest method in any given 
case. 



1.00 

10 
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(2) When k is large, generation of ‘-Kj exponential 
variates may require excessive time. In this case, taking 
the logarithm of the product of L k- f uniform variates may be 
a more efficient way to proceed. 

It may be possible to overcome both these 
difficulties by extending the decomposition technigue to 
values of k greater than one. In this case, however, the 
density function is no longer monotone and automatic 
generation of the sub-distributions and their probabilities 
will be considerably more difficult. 

2 . Tests for Tr en d 

As mentioned previously, a more thorough 

investigation of the theoretical properties of the Y test 

is also indicated. Its relative power compared to other 
tests against the trend model should be determined; possible 
examples of such tests include a normal theory regression 
test for trend after a logarithmic transformation or the 
various non- parametric tests for trend that are essentially 
equivalent to Kendall's rank correlation test. Results 
should also be obtained for other types of renewal process. 
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c 

cccccccccccccccccxccccccccccccccccccccccccccccccccccccccc.ccc 

SUBROUTINE GMI NIK K, BETA) 

ADDITIONAL ENTRY POINT: 

GAMA ( Z , IX) 

PURPOSE : 

GENERATION OF GAMMA RANDOM DEVIATES WITH SHAPE 
PARAMETER LESS THAN ONE. 

PROGRAMMER: 

D.W. ROBINSON, LT , USN 
METHOD : 

A MODIFICATION OF MARSAGLIA'S RECTANGLE -WEDGE-TA I L 
METHOD FOR NORMAL DEVIATES IS USED. THE GAMMA PDF 
IS DECOMPOSED INTO A HEAD REGION, A NUMBER 
(DEPENDENT ON K) OF RECTANGLES AND WEDGES AND A TAIL 
REGION. THE GMINIT SECTION CF THE SUBROUTINE ALSO 
SETS UP A BINARY SEARCH TREE TO BE USED FOR 
EFFICIENT SELECTION OF THE PFOPER REGION DURING THE 
ACTUAL GENERATING PROCESS, WHICH IS HANDLED BY THE 
GAMA SECTION. 

DESCRIPTION OF PARAMETERS: 

K GAMMA DISTIRBUTION SHAPE PARAMETER (MUST BE 

.GE. 0.05 AND .LT. 1.0) (REA L-4 ) 

BETA GAMMA DISTRIBUTION SCALE PARAMETER ( REAL* A ) 

IX SEED FOR UNIFORM RANDOM NUMBER GENERATOR 
Z RETURNED GAMMA DEVIATE (RtAL*4) 

NOTE THAT THE DENSITY FUNCTION OF A GAMMA RANDOM 
VARIABLE IS GIVEN BY: 

K K-L (-X/BETA) 

F ( X ) = (I/BETA) X E / GAMMA (K) 

THE FOLLOWING SUBROUTINES ARE CALLED: 

IGAM(KtX) COMPUTES THE INCOMPLETE GAMMA 

FUNCTION (GAMMA CDF). 

I N VGAM ( K , X ) COMPUTES THE INVERSE GAMA CDF 

RANDOM ( IX t U, N ) RETURNS A VECTOR U OF N UNIFORM 

RANDOM NUMBERS. 

OVFLOW SET-UP ENTRY POINT FOR RANCCM 

NOTE: 

UNDERFLOW IS POSSIBLE WHEN K IS LESS THAN .18 AND 
BECOMES MORE LIKELY AS K DECREASES. WHEN K IS .05 
THE PROBABILITY OF UNDERFLOW IS ABOUT .0001.3 FOR 
ANY GENERATED DEVIATE. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



c 



c 



SUBROUTINE GMINIT ( K , B ETA ) 

REAL*4 K, IGAM, INVGAM 
I NTEGE R*4 F I RST , TABLE » BOTTOM t END 
LOGICAL*! USED 



DIMENS I 
D I MENS I 
DIMENSJ 
DIMENSI 
ECU1VAL 



ON 


X( 101 


) , H 


( 1 


ON 


PROB ( 


202 


) r 


ON 


TEST ( 


202 


) , 


ON 


RAND ( 


2 ) 




ENCE (U, 


RAN 


D( 



00) ,P( 100) , Q( 100) , R( 
TA8LE( 20? ) , NEXT (202) 
LIST (302) , U S E D ( 202 ) 

1 ) ) , ( V , RAN D ( 2 ) ) 



100) ,B( 100) 
, LAST ( 202) 
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ThIS FIRST SECTION INITIALIZES CONSTANTS AND TABLES 
TO BE USED BY GAMA WHEN IT IS CALLED. THE FOLLOWING 
ARE TO BE CALCULATED: 

FIRST STARTING POINT FOR BINARY SEARCH 
NEXT, LINK VECTORS FOR BINARY SEARCH PROCEDURE 
LAST 

PROB VECTOR OF CUMULATIVE PROBABILITIES 
TABLE USED TO LOOK-UP RESULT OF SEARCH IN PROB 
J I POSITION IN PROB OF PO 

X(I) LEFT-HAND BOUNDARY OF i-TH RECTANGLE 
H ( I ) WIDTH OF I-TH RECTANGLE 
R ( I ) REJECTION RATIO FOR I-TH WEDGE 
B ( I ) Y ORDINATE FOR I-TH WEDGE 
ALPHA K— 1 

PO PROBABILITY FOR HEAD F EG I ON 

PN PROBABILITY FOR TAIL REGION 

HI TO CONSTANTS FOR HEAD APPROXIMATION TO INVERSE 
H4 GAMMA CDF 



CHECK FOR K IN RANGE 

IF((K .GE. 0.051 .AND. (K .LT. 1.0)) GO TO 5 
WRITE (6 ,4) K 

4 FORMAT! //' OGMINIT CALLED WITH K=',LPEL6.6, 

* ' OUT OF RANGE*/) 

RETURN 

GET UPPER BOUND ON NUMBER OF RECTANGLES REQUIRED 

5 N = 20. + 6.6/K 

I F ( N ,GT. LOO) N = LOO 

INITIALIZE CONSTANTS 

M = 2*N + 2 
MM = M - L 
ALPHA = K - L. 

GK - GAMMA(K) 

PO = 5.0E-5 / (K * K) 

HFAC = 2. 

FIND CONSTANTS FOR HEAD APPROXIMATION 

HI = K * GK 
H 2 = L < / K 

H3 = 2. 2E- 7 * (K + L. ) 

H 4 = 4 . / (K + L . ) 

SET UP RECTANGLE BOUNDS 

X ( L ) = (HL * PO) ** H2 

X ( L ) = 2. * X ( L ) / (L. + SQRT ( L . - H4 * X(LH) 

PO = I G AM ( K , X ( L ) ) 

HI = HL * PO 
H ( L ) - .2 5 * X ( L ) 

DO LO 1=2, N 

X ( I ) = X ( I - L ) + H( I-L) 

H ( I ) = H( I-L ) * FF AC 
P ( I ) = 0. 

Q ( I ) =0. 

LO CONTINUE 

X ( N+ L ) = X ( N ) + H ( N ) 

ZERO PROBABILITY VECTORS AND LINKS 

DO L5 I =L , M 

N E XT ( I ) = 0 
LAST { I ) = 0 
PROB(J) = 0 
LIST! I ) = 0 
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USED(I) = .FALSE. 

L 5 CONTINUE 

FIND COEFFICIENTS FOR NEWTON-R APFSON APPROXIMATION 
TO DENSITY FUNCTION TANGENT 



B I = 


-2. * 


ALPHA 


B 2 = 


ALPHA 


* ( ALPHA - L . ) 


AL = 


BL + 


1 . 


A 2 = 


ALPHA 


* ( ALPHA - 2 . ) 



C = 1. - ALPHA 

FIND RECTANGLE PROBABILITIES AND WEDGE VALUES 
PL = PO 

FL = E XP ( ALPHA * ALOG ( X ( 1 i I - X(L>) / GK 
DO AO I = 1 1 N 

FU = EX P ( ALPHA * AL0G(X(I + 1)) - X ( I + L ) ) / GK 
PU = I GAM ( K , X ( I + 1 ) ) 

P ( I ) = H< I ) * FU 
Q ( I ) = PU - PL - P( I ) 

NEWTON-RAPHSON ITERATION TO FIND POINT WHERE 
TANGENT TO PDF IS PARALLEL TO CHORD 

W = X ( I ) 

S = ( FU - FL ) / H ( I) 

SC = S * GK 
DO 20 J = L f 15 

Y = W * ( ( W + A L ) W + A2 + 

* SC * EX P ( C * ALCG(W) + W) ) 

* / ( ( W + B L ) * W + B 2 ) 

I F ( ABS ( Y - W) .LT. 1 . E-4 * Y> GO TO 30 
W = Y 

20 CONTINUE 

FIND VALUES FOR REJECTION METHOD 

30 B ( I ) = EXP ( AL PHA * ALOG(Y) - Y )/ GK + S*(X(I) - Y) 

R ( I ) = ( B ( I ) - FU) / (FL - FU) 

TEST TO SEE IF ENOUGH RECTANGLES HAVE BEEN TAKEN 
I F ( PU .GT. 0.999) GO TO 45 

RESET PROBABILITY AND FUNCTION VALUES FOR 
NEXT RECTANGLE 

PL = PU 
FL = FU 
40 CONTINUE 

FIND LOWER END INDEX FOR PR03 
45 LOW = 2 * (N-I ) + L 
FIND TAIL PROBABILITY 
PN = L. - PU 

GENERATE PROBABILITY VECTOR 

DO 80 1=1, N 

PROB ( I ) = P { I ) 

PROB ( I + N ) = G ( I ) 

TABLE ( I ) = I 
TA BL E ( I+N) = -I 
80 CONTINUE 

PROB(M-l) = PO 
PROB ( M ) = PN 
TABLEtM-l ) = 0 
TABLE(M) = 0 
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SORT PROBABILITY VECTOR 



0) .AND. ( PROB ( I ) .EC. P0 I 1 J 1 = 1 
+ PROB ( I-I ) 

5) .AND. (FIRST .EQ. 0))FIRST=I 



DC 110 1=1, MM 
ICH = 0 
L = M - I 
DO 100 J=1,L 

I F ( PRO B ( J ) - P R 0 8 ( J -J- 1 ) ) 100,100,90 

90 TEMP = PR0B(J) 

PR0B(J) = PROB(J-H) 

PROB ( J + 1 ) = TEMP 
ITEMP = TABLE(J) 

TABLE(J) = TABLEU + l) 

TABLEU+1) = ITEMP 
ICH = 1 

100 CONTINUE 

I F C ICH* 120,120,110 
110 CONTINUE 

CONVERT PROB TO CUMULATIVE PROBABILITIES 
FIND FI RST AND J 1 

120 Jl = 0 

FIRST - 0 
L = LOW + 1 
DO 130 I = L , M 

I F ( ( T ABL E ( I 5 .EQ. 

PROB(I) = P ROB ( I J 
I F ( ( P RO B ( I ( , GE . 

130 CONTINUE 

I F( FIRST .EQ. M) FIRST = MM 
PROB(M) = 1. 

I F ( J l .EQ. 01 Jl = LOW 

NOW DETERMINE THE VECTORS NEXT AND LAST 

150 BOTTOM = 1 
END = 1 

D D — 9C 

LIST! 1 ) = FIRST 

TEST(l) = .5 

USED ( F I RST ) = .TRUE. 

FIND NEXT ( K) AND L AST ( K ) FOR EACH K ON LIST 

1 5 L DO 15 9 1 = 1, BOTTOM 
LI = LI ST ( I > 

FIND NEXT { L I ) 

IF Ll+l HAS BEEN TESTED GO ON TO FINC LAST(LI) 

I F ( US ED ( LI + 1 ) ) GO TO 155 

GET PROBABILITY VALUE FOR NEXT ( L I ) 

PRN = TEST ( I ) + PR 

FIND J SUCH THAT PROB(J) > PRN 



152 

153 
l 5 A 



DO 152 J=LOW , MM 

I F ( PRO B ( J ) .GT. PRN) GO TO 153 
CONT INUE 

GET K SUCH THAT PROB(K) HAS NOT BEEN TESTED AND 
LI .LT. K ,LE. J. IF K EXISTS, SET NEXT(LI)=K. 

I F ( . NOT . US ED ( J ) ) GO TO 154 

J = J - 1 

I F { L I - J) 153, 155,155 
NEXT(LI) = J 
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155 



156 

157 

1571 

158 

1585 

159 

1591 

160 

161 

162 



NEXT (LI) HAS NOW BEEN FOUND IF IT IS NCN-ZERO. 

NOW ADD IT TO THE LIST IN PLACE OF LI FOR THE 
NEXT PASS THROUGH THE LIST. 

L I ST ( I ) = NEXT (L I ) 

NOW FIND LAST (LI ) 

IF LI = LOW CR PRQB(LI-l) HAS BEEN USED, GO 
ON TO THE NEXT VALUE ON THE LIST 

I F ( ( L I ,EQ. LOW) .OR. USE D ( L I— 1 ) ) GO TO 159 

GET PROBABILITY VALUE FOR LAST (LI) 

PRL = TEST ( I 1 - PR 

RESET PROBABILITY TEST VALUE FOR NEXT PASS 
TE ST ( I ) = PRN 

FIND J SUCH THAT PROB(J) > PRL 
DO 156 J = L0 W , MM 

I F ( P RO B ( J ) .GT. PRL) GO TO 157 
CONTINUE 

FIND K SUCH THAT PROB(K) HAS NOT BEEN TESTED AND 
LI .GT. K .GE. J. IF K EXISTS, SET LAST(LI)=K. 

IF ( .NOT .US ED ( J ) ) GO TO 157 1 
J = J + 1 

IFILI - J) 158, 158, 157 
LAST ( LI ) = J 
GO TO 1585 
J = LI 

LAST (LI) HAS NOW BEEN FOUND. ADD IT TC THE END CF 
THE LIST FOR THE NEXT PASS. 

END = END + 1 
L I ST ( END ) = J 
TEST(END) = PRL 
CONTINUE 

NOW RESET THE LIST FOR THE NEXT PASS BY ELIMINATING 
ZERO ENTRIES. 

BOTTOM = END' 

PR = PR * .5 
I = 1 

TEST FOR ZERO ENTRY 
IF(LIST(ID 160,160,163 

ZERO ENTRY FOUND. ADJUST BOTTOM OF LIST 

BOTTOM = BOTTOM - 1 

IF LIST EMPTY, QUIT. 

I F( BOTTOM) 165,165,161 

SHIFT OTHER LIST ENTRIES UP 

DO 162 J=I , BOTTOM 

L I ST ( J ) = L I S T t J+l) 

T E ST ( J ) = TEST (J + l ) 

CONTINUE 
GO TO 1591 
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C SET USED FLAG FOR LIST ENTRY 

C 

1.63 USED! L I ST ( I ) ) = .TRUE . 

GET NEXT LIST ENTRY 
1=1 + 1 

IF (I ,LE. BOTTOM) GO TO 1591 

DONE WITH LIST RESETTING. GO BACK FOR ANOTHER PASS 

END = BOTTOM 
GC TO 151 

SETUP FIRST CALL TO RANDOM 
165 CALL OVFLOW 

THE FOLLOWING STATEMENTS PRINT OUT A TABULAR LSTING 
OF THE RESULTS OF GMINIT. THEY MAY BE REMOVED, IF 
DESIRED, IN A PRODUCTION RUN. 

WRITE! 6, 170 IK, BETA 

170 FORMAT ( 1 1 GENE RAT ED VALUES FOR K= ' , 1 PE 14 . 6 , 

* ' BETA= ' , E 14. 6/ ) 

WRITE (6 , 175) PO, PN 

175 FORMAT ( ' OHEAD PROB AB I L I TY = • , 1 P E 14 .6 , 

* ' TAIL PROBABILITY^ ,E14.6//) 

WRITE ( 6 ,130) 

180 FORMAT! ’ OR EC T ANG L E / WEDGE V AL UE S ' / / 2X , • I • , 9X , ’ X ( I) 1 , 

4 12X, ' P( I ) • , 12X, *Q(I) r , L2X, * R ! I ) • , 12X, 1 B( I) ' /) 

ADD UP TOTAL RECTANGLE AND WEDGE PROBABILITIES AND 
PRINT OUT RESULTS 

SUM 1 = 0. 

SUM 2 = 0. 

DC 192 1=1, N 

NO ACTION FOR ZERO PROBABILITIES 

IF(PEI)) 193,193,185 

135 SUM1 = SUM 1 + P ( I ) 

SUM2 = SUM2 + Q ! I ) 

WRITE (6, 190) I ,X! I ) ,P( I ) ,Q( I ) ,R(1 ) , B( I) 

190 FORMAT! 1X13 , 1P5E16. 6) 

192 CONTINUE 

193 WRITE(6 ,194) SUM1 ,SUM2 

194 FORMAT! «OSUMS' , 1 5X , IP 2E 16 . 6 ) 

WRITE OUT CONSTANTS 

WRITE (6, 195) J1,H1,H2,H3,H4 

195 FORMAT! /• OVALUES FOR HEAD/TAIL A P PRCX I M AT I CN : » / 

4 8 J 1= 1 , 1 4/ • HL=',E16.6, ‘ H2='»E16.6, 

* * H3= ' , E 16 . 6 , ' H4= ' , E 1 6 . 6 ) 

WRITE OUT BINARY SEARCH DATA 
WRITE! 6 , 196) F IRST 

196 FORMAT ( / ' OSTART I NG POINT FOR BINARY SEARCH', 14) 

WRITE ( 6, 197) ( I , PROB ( I ) , TABLE! I ) , NEXT! I ) , LAST ( I ) , 

* I = L 0 W , M ) 

197 FORMAT (/• 0PR0BA6IL ITY VECTOR AND LINKS'// 

* 14 X , ' PROB ' , 9X, * TABLE * ,4X , * NEXT' , 5X , 

* 'LAST ' // ( IX I 15, 1PE 16.6, 0P3I 9) ) 

FINISHED WITH SET-UP PHASE. QUIT. 

RETURN 
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ThIS IS THE SECTION WHICH ACTUALLY GENERATES THE 
GAMMA VARIATES. 

ENTRY GAMA ( I X * Z ) 

GET TWO UNIFORM RANDOM VARIATES 

CALL RANDOM! IX, U, 2) 

CONDUCT BINARY SEARCH USING THE FIRST UNIFORM VARIATE 
J = FIRST 

200 I F ( U — PROB(J)) 210,250,230 

U < CURRENT VALUE. USE LAST FOR FOLLOWING TEST. 

210 I F ( LAST ( J ) ) 250,250,220 
220 J = LAST(J) 

GO TO 200 

U > CURRENT VALUE. USE NEXT FOR FOLLOWING TEST. 

230 I F ( NE XT ( J ) ) 245,245,240 
2 4 C J= NEXT(J) 

GO TO 200 

245 J= J + 1 

LOCATED PROPER PROBABILITY DIVISION. LOOK UP IN TABLE 

250 N = TABLE(J) 

I F < N I 260,290,320 

TABLE VALUE < 0. SAMPLE FROM N-TH WEDGE. 

260 N= -N 

GET ONE MORE UNIFORM DEVIATE 
CALL RANDOM! IX, U, I) 

GET U .LE. V 

270 I F ( U .LE. V) GO TO 280 
TEMP = U 
U = V 
V = TEMP 

GET TRIAL GAMMA VALUE 
280 Z = X ( N ) + H(N)*U 

PERFORM REJECTION TEST 
I F ( V .LE. R ( N I ) GO TO 330 

FIRST TEST FAILED. GET VALUE FOR SECOND TEST. 

W = U + EXP! ALPHA * ALOG(Z) - Z) / B ( N ) 

I F ( V . L E . W ) GO TO 330 

SECOND TEST ALSO FAILED. REPEAT LOOP WITH TWO 
NEW UNI FORM VARIATES 

CALL RANDOM! IX, U, 2 ) 

GC TO 270 

TABLE VALUE = 0. USE HEAD/TAIL DISTRIBUTIONS. 
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C DETERMINE WHICH DISTRIBUTION TO USE. 

C 

290 I F ( J .EQ. Jl) GO TO 300 

USE TAIL DISTRIBUTION. 

Z = INVGAM (Kf PN * V) 

GC TO 330 

USE HEAD DISTRIBUTION. 

300 Z = (HI * V) ** H 2 

I F ( Z .LT. H3 ) GO TO 330 
Z = 2. * Z / (L. + SQRT ( L . - H4 * Z > ) 

GO TO 330 



TABLE VALUE > 0. SAMPLE FROfr N-TH RECTANGLE. 
320 Z = X( N ) + H ( N I * V 



Z IS NOW A UNIT GAMMA VAR I AT I: . SCALE IT AND EXIT. 

330 Z = Z * BETA 
RETURN 
END 
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GAMA 



CSECT 





ENTRY 


START 




BASE 1 


EQU 


2 




BASE2 


EOU 






SAVE 


EQU 


4 




REGA 


EQU 


5 




REG B 


EQU 


6 




REGC 


EQU 


8 




SRC H 


EQU 


7 




£ 


STM 


14, REGC , 1 2 ( 13) 


SAVE CALLING PROGRAM REGS 




BALR 


BASEL , 0 


SET UP BASE REGISTERS 




USING 


*,BASE1 






LA 


BASE2, DATA 






USING 


DATA, BASE2 






LR 


SAVE , 13 


SAVE ADDRESS OF CALLING 




LA 


13 , SV AR EA 


PROGRAM SAVE AREA 
LOAD OWN SAVE AREA ADDRESS 




ST 


13 ,8( 0, SAVE) 


STORE OWN SAVE AREA ADDRES 




ST 


S AV E , 4 ( 0 , 13) 


STORE CALLING SAVE AREA 


* 


LA 


1, ARG 1 


ADDRESS 

LINK TO RANDOM FOR TWO 




L 


15, =V ( RANDOM ) 


UNIFORM DEVIATES 




BALR 


14,15 






L 


REGD, F 1RST 


START BINARY SEARCH 




L 


SRCH , U 


LOAD UNIFORM FCR SEARCH 


❖ 

SEARCH 


BALR 


REGC , 0 


LOAD BRANCH ADDRESS 


LR 


RE CA, REGB 


LOAD NEXT INDEX 




C 


SRCH, PROB (REGA) 


COMPARE 




BC 


4 , L OW 






BC 


8, FOUND 




HIGH 


L 


RFGB, NEXT( REGA) 


GET NEXT INDEX FROM "NEXT" 




BCTR 


REGB, REGC 


BRANCH BACK IF INDEX 




LA 


REGA, 4( REGA) 


IS NON-ZERO 


a. 


BC 


L5, FOUND 


THROUGH WITH SEARCH 


LOW 


L 


REGB, LAST (REGA) 


GET NEXT INDEX FROM "LAST" 




BCTR 


REGB, REGC 


BRANCH BACK IF NON-ZERO 


FOUND 


L 


REGB, TABLE ( REGA ) 


FIND WHICH METHOD TO USE 




LTR 


REGB, REGB 






BC 


2, RECT 


RECTANGLE IF TABLE 


JU 

a> 






VALUE GREATER THAN 0 


*v> 


STD 


2 , SAVED 


SAVE FLOATING POINT REG 2 




BC 


8 , TAI L 


USE HEAD/TAIL IF TABLE 


❖ 






VALUE = 0 


WEDGE 


LCR 


REGB, REGB 


USE WEDGE METHOD IF TABLE 


❖ 




• 


VALUE < 0 




LA 


1 , ARG 2 


GET NEW RANDOM NUMBER 




L 


15, =V (RANDOM) 




* 


BALR 


14,15 




WL 


LE 


0,U 


GET U LESS THAN OR EQUAL V 




CE 


0, V 




* 


BC 


L 2 , W 2 




LE 


2 , V 


EXCHANGE U AND V 




STE 


0, V 






LER 


0,2 




* 


STE 


0,U 




W2 


ME 


0 , H ( R EG B ) 


FIND Z = X ( N ) + U * H ( N ) 




AE 


0 , X ( R EG B ) 
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LE 


2, V 




CE 


2, R( REGB) 




BC 


12 , DONE 




STE 


0»Z 




LA 


1 , ARG3 




L 


15,=V( ALOG! 




BALR 


14 , 15 


* 




ME 


O, ALPHA 




SE 


o,z 




STE 


0, ZP 


* 




LA 


1 , ARG3A 




L 


1 5 , = V ( E XP ) 




BALR 


14, 15 


* 




DE 


0 , B ( REGB) 




AE 


0,U 




CE 


0, V 




LE 


0,Z 




BC 


2 , DONE 


❖ 




LA 


1 , ARG 1 




L 


15, =V (RANDOM) 




BALR 


14,15 




BC 


1 5 , W 1 


* 


TAIL 


C 


REGA, J1 




BC 


8, HEAD 




LE 


0, PN 




ME 


0, V 




STE 


o,z 




LA 


1 , ARG4 




L 


1 5 , = V ( I NVGAM) 




BALR 


14, L5 




BC 


15, DONE 


* 


HEAD 


LE 


0, V 




ME 


0 , HI 




STE 


o,z 




LA 


1 , ARG5 




L 


15,=V( FRXPR# ) 




BALR 


14, 15 




CE 


0, H3 




BC 


4, DONE 




STE 


OrZ 




ME 


0 , H4 




LE 


2 » ONE F L 




SER 


2,0 




STE 


2 ,H5 




LA 


1 , ARG6 




L 


L 5 , = V ( SQRT) 




BALR 


14,15 




AE 


0,0NEFL 




LE 


2 , Z 




ME 


2 , TWOFL 




DER 


2,0 


* 


LER 


0,2 


DONE 


LD 


2, SAVED 




BC 


15, SCALE 


* 


RECT 


LE 


0 , H( R EG B ) 




ME 


0, V 




AE 


0,X( REGB) 


* 


SCALE 


ME 


0, BETA 




LR 


13, SAVE 


<A> 



TEST IF V IS LESS THAN OR 
EQUAL TO R(N) 

PASS REJECTION TEST IF SO 
THIS STEP IS REACHED ONLY 
1/3 OF THE TIRE 
FIND ALOG(Z) 



ZP = (K - 1 ) * ALOG ( Z ) - Z 

FIND EXP(ZP) = 

EXP(-Z) * (Z ** K-l) 



FIND W = U + 

EXP(-Z)-( Z ** K-l ) /B (N) 
TEST IE V IS LESS THAN OR 
EQUAL TO U 
QUIT IF SO 

OTHERWISE GET TWO NEW 
UNIFORM DEVIATES 

AND REPEAT THE LOOP 

SEE IF HEAD OR TAIL IS 
REQUIRED 

TAIL WILL BE USED. 

Z = I NVGAM ( K , PN*V) 



COMPUTE FIRST ORDER APPROX 
( V * PO * GAMMA ( K+ 1 ) ) 

* # 1 / K 



SEE IF SECOND ORDER APPROX 
IS NECESSARY 

Z = 2*Z / (1 + 

SQRT ( 1 - 4*Z/ (K+l ) ) ) 



RESTORE FLOATING POINT 
REGISTER 2 



Z = X ( N ) + V * H ( N ) 



SCALE VARIATE AND EXIT 

RESTORE OLD SAVE AREA 
ADDRESS 
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LM 


2 i R EGC j 28 ( 13) 


RESTORE OLD GENERAL 


* 


L 


14,12(13) 


PURPOSE REGISTERS 
PESTCRE RETURN ADDRESS 




MV I 


12(13) tX'FF* 


SET END OF SUBROUTINE FLAG 


a# 


BCR 


15, 14 




V 


DROP 


BASE1 




£ 


DROP 


BASE2 




Itart 


ENTRY POINT FOR INITIALIZING ROUTINE 


STM 


14 ,REGC , 1 2 ( 13) 


SAVE GP REG I STERS 




BALR 


BAS El , 0 


SET UP BASE REGISTERS 




USING 


* , BA S E 1 






LA 


BASE2, DATA 






USING 


DATA, BASE2 






LR 


SAVE, 13 






LA 


13, SVAREA 






ST 


13,8(0, SAVE ) 






ST 


SAVE,4( 0, 13) 






L 


REGA, 0( 1 ) 


GET SET UP ARGUMENTS 




MVC 


K ( 4 ) , 0( REGA) 


K (SHAPE PARAMETER) 




L 


REGA, 4( L) 






MVC 


B E T A ( 4 ) , 0 ( REGA ) 


BETA (SCALE PARAMETER) 


A 


L 

MVC 


REGA, 8( 1) 

I X ( 4 j , 0 (REGA) 


IX (RANDOM NUMBER SEED) 


T 


LA 


1 , ARG7 


CALL GMINIT 




L 


15,-V(GMINIT) 




A 


BALR 


14,15 




*T 


L 


15,=V(0VFL0W) 


CALL OVFLOW TO SET UP FOR 


* 


BALR 


14,15 


RANDOM 




LR 


13, SAVE 






LM 


2 , REGC ,28(13) 






L 


14, L 2 ( 13) 






MV I 


12(13) ,X'FF‘ 




* 


BCR 


15,14 




* 


DS 


OD 




IX 


DC 


4X * 00 ' 




u 


DC 


4 X s 00 ' 




V 


DC 


4X e 00' 


*■ 


z 


DC 


4X e 00 ' 




ZP 


DC 


4X* 00 1 




ONE 


DC 


F* 1 ' 




TWO 


DC 


F ' 2 * 




ONE F L 


DC 


E ' 1 .0 ' 




TWOFL 

A 


DC 


E * 2 • 0 ' 




"v 

K 


DC 


4X 1 00 1 




BETA 

A 


DC 


4X' 00' 




AL PHA 


DC 


4 X * 0 0 1 




FIRST 


DC 


4X* 00' 




PR06 


DC 


4 X ' 00* 




PI 


DC 


808X' 00 • 




NEXT 


DC 


4X' 00 s 




NL 


DC 


808X S 00 ' 




LAST 


DC 


4X» 00 • 




LI 


DC 


8 08 X ' 00' 




TABLE 


DC 


4X‘ 00 ' 




T 1 


DC 


808X' 00* 




DATA 


DS 


OF 




X 


DC 


4X' 00* 




XI 


DC 


404X' 00 • 




H 


DC 


4X' 00 • 




HP 


DC 


4 0 0 X * 00' 




R 


DC 


4X e 00 • 
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R 1 


DC 


40 OX ' 00 * 


B 


DC 


4 X ' 0 0 1 


B 1 


DC 


400X' 00 • 


JL 


DC 


4 X * 00 ' 


PN 


DC 


4X* 00* 


HI 


DC 


4 X * 00 ' 


H2 


DC 


4X* 00 • 


H3 


DC 


4X ' 00 * 


H4 


DC 


4X ' 00 * 


H5 


DC 


4 X ' 00 • 


argl 


DC 


A( IX) 




DC 


A ( U ) 




DC 


X * 80 ' 




DC 


AL3 ( T WQ ) 


ARG2 


DC 


A( IX) 




DC 


A( U) 




DC 


X * 80* 




DC 


AL3 ( ONE ! 


ARG3 


DC 


X’ 80* 




DC 


AL3 ( Z ) 


ARG3A 


DC 


X* 80' 




DC 


AL3( ZP ) 


ARG4 


DC 


A ( K ) 




DC 


X' 80’ 




DC 


AL 3 ( Z ) 


ARG5 


DC 


A ( Z ) 




DC 


X' 80* 




DC 


A L 3 ( H 2 ) 


ARG6 


DC 


X * 8 0 * 




DC 


A 13 ( H 5 ) 


ARG7 


DC 


A (K ) 




DC 


A ( ALPHA ) 




DC 


A ( FIRST ) 




DC 


A ( P L ) 




DC 


A ( N I ) 




DC 


A ( L 1 ) 




DC 


A ( T I ) 




DC 


A ( X I ) 




DC 


A ( HP ) 




DC 


A ( R L ) 




DC 


A ( B I ) 




DC 


A( Jl) 




DC 


A( PN ) 




DC 


A ( H 1 ) 




DC 


A ( H2 ) 




DC 


A ( H3 ) 




DC 


X * 80 * 




DC 


AL3 ( H4 ) 


SVAREA 


DS 


L8F 


SAVED 


DS 


3D 




END 





on o on ooooo on on on no no on no on on ooooooooo ooooooo 



FUNCTION I GAM ( K , X ) 

PURPOSE : 

EVALUATION OF THE INCOMPLETE GAMMA FUNCTION: 

K-l -Y 

INTEGRAL OF Y E FROM 0 TO X DIVIDED BY 
GAMMA! K) 

METHOD : 

FOR VALUES OF X LESS THAN 5 AN INFINITE SERIES 
APPROXIMATION IS USED; FOR X GREATER THAN 5 A 
CONTINUED FRACTION APPROXIMATION IS USED. IF X 
IS LESS THAN OR EQUAL TO ZERO, ZERO IS RETURNED. 



IMPLICIT REAL*8(D) 

REAL-4 I GAM , K , V ( 6 ) 

REAL-8 EPS/I. D-I3/ 

TEST VALUE OF X 
I F ( X .GT. 0 . ) GO TO 10 

X .LE. 0. RETURN IGAM = 0. 

I GAM= 0 . 

RETURN 

DECIDE WHICH APPROXIMATION TO USE 
10 I F ( X .GT. 5 . 0 ) GO TO SO 

USE INFINITE SERIES APPROXIMATION 
DX=DBLE (X) 

DK=DBLE (K ) 

INITIALIZE SUM AND TERM 
DTERM=DX**DK/DGAMMA(DK ) 

D SUM=DT ERM/DK 

SUM THE SERIES 
DO 30 1=1,30 

TEST FOR CONVERGENCE 
I Ft DABS(DTERM) -EPS*DSUM)40 ,40 ,20 

HASN'T CONVERGED. ADD ANOTHER TERM 
20 DTERM=DTERM*(-DX»/DFLOAT( I) 

DSUM=DS UM + DT.ERM/ (DK+D FLOAT! I ) ) 

30 CONTINUE 

SERIES CONVERGED. EXIT. 

40 I GAM=SNGL ( DSUM ) 

RETURN 



USE CONTINUED FRACTION APPROXIMATION 

INITIALIZE CONTINUED FRACTION VALUES 
50 Y= 1 , - K 
W=X+Y+L .0 

SET UP RATIO VECTOR 

V ( l ) = 1 . 0 
V ( 2 ) =X 

V! 3) =X+ 1. 0 

V (4) = W*X 

R = V ( 3 ) / V ! 4 ) 

COUNT = l .0 

GET NEXT APPROXIMATION 
60 Y= Y+ 1 . 0 
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W=W+2.0 

V (5 )=V ( 3 )*W-V( 1 ) *Y*COUNT 
V ( 6 ) = V ( 4)*W-V( 2)*Y*C0UNT 
R AT 1 0= V ( 5 ) / V ( 6 ) 

TEST FOR CONVERGENCE 

I F ( ABS ( RAT IO-R ) .LE. L.0E-6*R)G0 TO 90 

DIDN'T CONVERGENCE. SET UP FOR NEXT ITERATION. 
R = RAT I 0 

CCUNT=CCUNT+1 .0 
DO 70 1=1,4 
V ( I >=V ( 1 + 2 ) 

70 CONTINUE 

TEST IF SCALING REQUIRED 
I F ( V ( 5 ) ,LT. I.OE 35 ) GO TO 60 

SCALE DOWN RATIO VECTOR 
DO 80 1=1,4 
V ( I ) = V ( I ) *L.0E-35 
80 CONTINUE 
GO TO 60 

CONVERGED. CONVERT RESULT AND EXIT. 

90 IGAM = L.-(RATIO*EXP(K*ALOG(X)-X)/GAMMA(KM 
RETURN 
END 
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FUNCTION INVGAM ( K * Z ) 



THIS FUNCTION SUBPROGRAM SOLVES THE EQUATION 
G ( K ; X J = L - Z 

BY NEWTON-RAPHSON ITERATION, WHERE G £ K ; X ) IS THE 
INCOMPLETE GAMMA FUNCTION WITH SHAPE PARAMETER K LESS 
THAN L.O. 

REAL*4 K» INVGAM, KP,KBAR 
DATA EPS/L.OE-6/ 

DIMENSION V ( 6 ) 

TEST FOR Z IN RANGE 
IF (Z .GT. 0. ) GO TO 10 
I N V GA M = 0 « 

RETURN 

TEST IF NEW K VALUE 
10 I F ( K .EQ. KP) GO TO 20 

NEW VALUE. RECOMPUTE K-DEPENDENT CONSTANTS. 

K P = K 

K B A R= l . — K 
GK=GAMMA { K ) 

INITIALIZE ITERATION VALUES 
20 ZBAR=-Z*GK 
X=I .0 

MAIN ITERATION LOOP. 

INITIALIZE CONTINUED FRACTION VALUES. 

30 Y=KBAR 
W=X+Y+ 1 .0 
V ( 1) = 1 . 0 
V ( 2 > = X 
V (3) =X+i .0 
V( 4}=W*X 
R = V ( 3 ) / V ( 4 ) 

COUNT = I .0 

CONTINUED FRACTION APPROXIMATION FOR G £ K ; X ) 

40 Y = Y + L o 0 
W=W+2. 0 

V(5)=V( 3)*W-V( U *Y*C0UNT 
V(6) =V( 4)*W-V(2 )*Y*C0UNT 
RATI0=V(5)/V(6) 

TEST FOR CONVERGENCE 

I F ( ABS ( RATIO-R) ,LE. EPS^R) C-0 TO 70 

DIDN'T CONVERGE. SET UP FOR NEXT ITERATION. 

R = RAT I 0 

CCUNT=COUNT+L .0 
DO 50 I=L,4 
V ( I ) = V ( 1 + 2) 

50 CONTINUE 

TEST IF SCALING REQUIRED 
I F ( V ( 5 ) .LT. 1.0E 35 ) GO TO 40 

SCALE DOWN V VECTOR TO PREVENT OVERFLOW. 

DC 60 I =L ,4 
V( I )=V( I)*L.0E-35 
60 CONTINUE 
GO TO 40 

APPLY NEWTON'S FORMULA 
7 0 XNEW=X*(1.+RATI0+ZBAR*EXP(X-KP*AL0G{X) ) ) 

TEST FOR CONVERGENCE 

I F ( ABS ( XNEW-X ) .LE. EPS*XNEW) GO TO 80 
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X=XNE W 
GC TO 30 



NEWTON'S METHOD CONVERGED. SET FUNCTION VALUE 
AND EXIT. 

80 I NVGAM= XNE W 
RETURN 
END 
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INVGAM 


CSECT 




SAVE 


EQU 


2 


BASE 


EQU 


3 


REGA 


EQU 


4 


REG B 


EQU 


5 


REGC 


EQU 


6 




STM 


1.4, REGC, 1 2 ( L 3 ) 




LR 


BASE, 15 


u- 


USING 


I NVGAM, BASE 


-f 


LR 


SAVE, 13 




LA 


13, SVAREA 


* 




ST 


13 ,8(0, SAVE) 




ST 


SAVE, 4(0, 13) 




STD 


2, SAVED 




STD 


4, SAVED-:- 8 


o> 


STD 


6, SAVED+ 16 




L 


R EGB , 4 ( 1 ) 




LE 


0 , 0 ( R EG B ) 




LTER 


0,0 




BC 


2 , T 1 




LE 


0 , = E ' 0 . 0 ' 




LR 


13, SAVE 


❖ 




LM 


2, REGC, 28(13) 




L 


14, 1 2 ( 13) 




MV I 


12(13), X'FF* 




BCR 


14,15 


* 


T 1 


STE 


0, Z 


* 




L 


REGB , 0 ( 1 ) 




LE 


0 , 0 ( REG B ) 




CE 


0, K 




BC 


8 , T 2 




STE 


0 ,K 




LCER 


2,0 




STE 


2, KNEG 




AE 


2 , - E 1 1 . 0 1 




STE 


2, KBAR 




LA 


1 » ARG 1 




L 


15, =V< GAMMA i 




BALR 


14, 15 




STE 


0 , G K 


❖ 


T 2 


LE 


0,Z 




LCER 


0,0 




ME 


0, GK 




STE 


0, ZBAR 




LE 


0 , =E ' 1.0' 




STE 


0,X 


* 


MAIN LOOP FOR NEWTON 


* 


LOOP 1 


LE 


0, KBAR 




STE 


0 , Y 




LE 


2 , = E *1.0' 




STE 


2 , V 




LE 


4 , X 




STE 


4, V+4 




LER 


6,4 




AER 


4,2 




STE 


4, V+-8 




AER 


0,4 




STE 


0,W 




MER 


0,6 



REGISTER EQUATES 



SAVE GENERAL PURPCSE REGS 
SET UP BASE REGISTER 



SAVE CALLING PRCGRAM SAVE 
AREA ADDRESS 

PUT OWN SAVE AREA ADDRESS 
IN REGISTER 13 
STORE SAVE AREA ADDRESSES 

SAVE FLOATING POINT REGS 



GET ADDRESS OF 1 - Z 

LOAD I - Z 

TEST IF I - Z > 0.0 

Z .LE. 0; SET INVGAM TO 0 
AND RETURN 

RESTORE CALLING PROGRAM 
SAVE ADDRESS 
RESTORE GP REGISTERS 
RESTORE RETURN ADDRESS 
SET END OF SUBROUTINE FLAG 
RETURN 

SAVE VALUE OF l - Z 

GET ADDRESS OF K 
LOAD VALUE OF K 
TEST IF VALUE IS NEW 

VALUE IS NEW. COMPUTE NEW 
CONSTANTS . 

KNEG = -K 

KBAR = 1 - K 

GET GAMMA ( K ) 



GK = GAMMA! K) 

COMPUTE Z-RELATED 
CONSTANTS 

ZBAR = -GAMMA(K) * (I - Z ) 
X = I 

RAPHSCN ITERATION 

INITIALIZE CONTINUED 
FRACTION CONSTANTS 
Y = 1 - K 

V( l) = 1.0 

V( 2) = X 



V ( 3 ) = X + 1.0 
W = X + Y + 1 . 0 
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STE 0 » V+ 1 2 V ( 4 ) = W * X 

DER 4,0 

STE 4 , PROB PROB = V(3) / V(4) 

STE 2 » NUM HUM = L 

CONTINUED FRACTION APPROXIMATION 



* 


L00P2 


LE 


0 , Y 




AE 


0 , =E * L.O* 




STE 


0,Y 




LE 


4, W 




AE 


4 , =E ' 2 . 0 ' 




STE 


4, W 




MER 


2,0 


* 




LE 


0, V+12 




MER 


0,4 




LE 


6, V+4 




MER 


6,2 




SER 


0,6 




STE 


0, V+20 




LE 


0, V+8 




MER 


0,4 




LE 


6, V 




MER 


6,2 




SER 


0,6 




STE 


0 , V+ l 6 


❖ 




LER 


4,0 


* 




DE 


0, V+20 




LER 


2,0 




SE 


0, PROB 




STE 


2, PROB 




LPER 


0,0 




ME 


2 » = E * I .OE-6' 




CER 


0,2 




BC 


12 , CNVRG 




MVC 


V( 16) , V+8 




CE 


4 , = E ' l . OE 3 5 ‘ 




BC 


12, INCMNT 


* 




LE 


0 , = E * l .OE-35 




LER 


2,0 




ME 


0, V 




STE 


0,V 




LER 


0,2 




ME 


0, V+4 




STE 


0, V+4 




LER 


0,2 




ME 


0, V + 8 




STE 


0, V+0 




LER 


0,2 




ME 


0, V+ 12 




STE 


0, V+12 


* 


INCMNT 


LE 


2, NUM 




AE 


2 , = E ‘ 1.0* 




STE 


2, NUM 




BC 


15, L00P2 


* 

* 


CNVRG 


LA 


l, ARG2 




L 


15, =V( ALOG) 




BA L R 


14, 15 






ME 


0 , KNEG 




AE 


0, X 


sjc 


STE 


0,XP 




LA 


1, ARG3 



Y = Y + L.O 



W = to + 2.0 



7(6) = V ( 4 ) * to - 

V( 2) * Y * NUM 



V ( 5 ) = V ( 3 ) * W - 

V ( 1 > * Y * NUM 

HEW RATIO = V ( 5 5 / V ( 6 i 
TEST FOR CONVERGENCE 

A B S ( P A T I 0 - PROS) : 

EPS * PROB 
SHIFT RATIO VECTOR 
TEST IF SCALING NEEDED 



SCALE DOWN RATIO VECTOR TO 
PREVENT OVERFLOW 



INCREMENT COUNT 



FIND ALOG(X) 



XP = -K * ALOG(X) + X 
FIND EXP(XP) 
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L 


1 5 , = V ( E XP ) 


* 


BALR 


LA, 15 


* 


ME 


0 , ZBAR 




AE 


0, PROB 


1 


AE 


0,=E • 1. O' 




ME 


0,X 


* 


LER 


2,0 




SE 


2 , X 




LPER 


2,2 




STE 


0, X 




ME 


0 , = E 1 L . OE-6 ‘ 




CER 


2,0 




BC 


2, LOOP! 


$ 


LE 


0,X 




LD 


2, SAVED 




LD 


4 , SAVED-!- 8 




LD 


6, SAVED+ L6 




LR 


13, SAVE 




LM 


2 , REGC , 28 ( 13) 




L 


L4, L2 ( L 3 ) 




MV I 


12(1. 3)jX c FF* 




BCR 


L 5 , L 4 


DATA 

* 


DS 


OD 


SAVED 


DS 


3D 


SVAREA 


DS 


L8F 


•V 

K 


DC 


4X * 00 1 


Z 

* 


DC 


4X * 00 » 


KNEG 


DC 


4 X ' 0 0 ' 


KBAR 


DC 


4 X ' 0 0 * 


GK 


DC 


4X' 00 * 


ZBAR 

* 


DC 


4X' 00' 


X 


DC 


4X' 00 ' 


V 


DC 


4X* 00 • 


V 


DC 


24X 1 00 * 


W 


DC 


4 X f 00 ‘ 


PROB 


DC 


4X' 00' 


NUM 


DC 


4X ' 00 e 


XP 

* 


DC 


4X* 00 ' 


ARG L 


DC 


X' 80 f 


$ 


DC 


AL3 ( K ) 


ARG2 


DC 


X ‘ 80 * 


4 


DC 


AL3 ( X ) 


ARG3 


DC 


X ' 80 ' 




DC 

END 


AL3 (XP) 



XNEW = X * (1 + RATIO - 

X — K 

( I-Z )*GAMMA(K) E X 
(BASIC NEWTON-RAPHSON 
RECURSION RELATION) 



TEST FOR NEWTON-RAPHSON 
CONVERGENCE 



ABS (XNEW - X ) : 

EPS * XNEW 

DONE. LOAD FUNCTION VALUE. 

RESTORE FLOATING POINT 
REGISTERS 

RESET OLD SAVE AREA 
RESTORE OLD GP REGISTERS 
RESTORE RETURN ADDRESS 
SET END OF SUBROUTINE FLAG 
RETURN 
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on oo oo on no no no on on 



* 



REAL*4 K,M2,MD 
REAL* 8 M,R,R1,DBLE 



R E AL*8 YL 
DIMENSION 
DIMENSION 
D IMENS ION 
DIMENSION 
D I MENS I ON 
CAT A 
DATA 



DATA 



»Y2 

M ( 4 ) , M2 (4) ,MD(4, LO ) 

K ( 4 ) 

QUANT (23) ,02(23) »QD( 23, 10) » SCORE (23) 

AL PH A ( 23 ) , I Q (23 ) 

S( 10000) ,U( LOO) 

NS AMP L/ 10000/ , SAMPL/ 10000./ 

ALPHA/ .001 , .002 , .005 , . 0 1. , .0 2 , .0 2 5 , .05 , 
. 2 i . 3 , .4? . 5 j . 6 , *7 , .8? .9 , .93, ©975 
.99, .995$ .998, .999/ 

K/.l, .25, .50, .75/ 



1 » 

. 98, 



INITIALIZE CONSTANTS 
I X= 135721 83 
DC 10 1=1,23 

IC(I ) =SAMPL* ALPHA ( I )+0 .5 
10 CONTINUE 
N = 30 



START SIMULATION LOOP 
DO L60 INDEX=l,4 

INITIALIZE GAMMA GENERATOR 
CALL START (K { INDEX ), 1 IX) 

FIND THEORETICAL MEAN AND VARIANCE 
TMEAN=( FLOAT (N)*L . )*.5 

TVAR= (F LOAT ( N) **2-1 . ) / ( 12 . * ( K ( I NDEX ) *Fl 0 AT ( N ) + 1 . ) ) 
T SD=SQRT ( T VAR ) 

ZERO OUT DATA COLLECTION ARRAYS 
DG 20 1=1,4 
M ( I ) =0 . 

M2( I)=0„ 

20 CONTINUE 

DC 30 1=1,23 
QUANT( I ) =0 . 

30 CONTINUE 

PERFORM 10 REPLICATIONS OF THE SIMULATION 
DO 90 L = L, 10 

GENERATE 10,000 RATIO STATISTICS 
DO 60 J = 1 , NS A MPL 

Y 1 = 0. 

Y 2 = 0 . 

DG 40 I =1 , N 
Y1=Y1 + DBLE(GAMA( IX ) > 

Y 2 = Y 2 -i- Y 1 
40 CONTINUE 

FIND RATIO 
R=Y2/Y1 
S ( J ) = R 
R L = R 

ACCUMULATE MOMENT ESTIMATORS 
DC 50 1=1,4 
M ( I ) = M ( I HR1 
R1=R1*R 
50 CONTINUE 
60 CONTINUE 

C DETERMINE MOMENT ESTIMATES FOR CURRENT SAMPLE 

DC 70 1=1,4 
MD { I , L ) =M( I ) /SAMPL 
M2( I ) = M2( I )+MD( I , L ) 

M ( I )=0. 

70 CONTINUE 
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c 



c 

c 



c 

c 



c 

c 



c 

c 



c 

c 



SORT SAMPLE AND GET QUANTILE ESTIMATES 
CALL S ORT ( S » 1,NSAMPL) 

DO 80 1=1.23 
QD( I,LI = S( IQ ( I )) 

QUANT ( I ) = QUAN1 ( I ) +S( I Q ( I ) ) 

80 CONTINUE 
90 CONTINUE 

FIND MEAN AND VARIANCE OF THE 10 MOMENT ESTIMATES 
DO 110 1=1,4 
M ( I ) = M2 (I )/lO . 

M 2 ( I) = 0 . 

DO 100 J= 1, 10 

M2 ( I) =M2 ( I ) + ( M( I ) - M D ( I, J) 1**2 
IOC CONTINUE 

M2 ( I J=M2< I J/9. 

110 CONTINUE 

FIND MEAN AND VARIANCE OF THE 10 QUANTILE ESTIMATES 

DO 130 1=1,23 

QUANT ( I )=QUANT( I ) / 10. 

Q2( I 1=0 . 

DC 120 J=1,10 

Q2( I J=Q2( I )+( QUANT ( I ) -QD( I , J ) J **2 
120 CONTINUE 

Q2 { I ) =02(1 ) / 9 « 

SCORE (I)=< QUANT ( I ) -TMEAN) /T SD 
130 CONTINUE 

COMPUTE SAMPLE STATISTICS 

Y ME AN= M ( 1 I 

V AR= SANPL*( M(2) -M( 1)*M(U)/(SAMPL-1.) 

SD=SQRT (VAR) 

SKEW= ( M(3)-( 2.*VAR+M( 2 ) )*M ( 1 ) )/ ( SC*VAR) 

C KURT= ( ( ( -3.*M(1) *M( l )+6.*M( 2) ) * M ( 1 )-4 . *M ( 3 ) ) 

* *M(l)+M(4) )/(VAR*VARJ 

PRINT OUT RESULTS 

WRITE (6, 140)K ( INDEX) ,N, M, M2 , YMEAN . VAR , SD , SKEW,CKURT 
140 FORMAT ( ‘ 1 SI MULAT ION RESULTS FOR K=',F5.2,‘ AND SAMPL', 

* *E SIZE = * , 14 

* ///'OTHE FIRST FOUR MOMENTS ARE : • / / L 0 X 1 P4E 1 9 . 6/ 

* 1 VARIANCE , 4E 19. 6/ // * 0 MEAN = ’, 014.6/ 

* ‘OVARIANCE =', F 14.6, 4X, ' STANDARD DEVIATIONS, 

* E14.6/ C 0SKEWNESS = ' , E 14 . 6/ 8 OKURTO S I S =*,E14.6/) 



WRITE (6 ,150). ( ALPHA { I) , QUANT ( I ) , Q2 ( I ) , SCORE ( I ) , I = 1 , 2 3 ) 
150 FORMAT { // / ' OQUANT I LES OF THE DISTRIBUTION:' 

* //5X , “ALPHA * , 6X, ' AL PH A- QUANT I L E * , 8X , ' VAR I ANCE ' 

* ,3X, 'NORMALIZED QUANTILE*// 

* (3X0PF7.3, 1P2E19.6,0PF14.6) ) 

160 CONTINUE 

STOP 

END 
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S L E ROUT IN 
D I NEKS I CN 
INTEGER’^ 
F = 1 



E SCRT ( A » I I » J J ) 

A ( J J ) t IL( 16) » I L ( 16) 
A » T , T T 



(A 

1C 



20 



2 C 



4C 

5C 



6 C 



1= I I 
J = J J 

I F ( I .GE. J) 
K= I 

I J= ( I + J )/2 
T= A ( I J ) 

IF(A{ I ) .LE. 
MI J)=A ( I ) 
A(I )=T 
T = A ( I J ) 

L = J 

IF(A(J) .GE. 
A ( IJ) = A(J) 

A ( J ) =T 
T = A ( I J ) 

I F ( A( I) .LE. 
A(IJ)=A(II 
A ( I ) =T 
T = A ( I J ) 

GC TO 4 C 
A ( L ) = A ( K ) 

A ( K ) = TT 
L = L-1 

I F ( A ( L i .GT. 
TT=A ( L ) 

K = K+l 

I F ( A ( K ) ,LT. 
I F ( K . L E . Li 
I F ( L- I .LE. 

T I ( W ) = T 

IL(F.Ul 
I = K 
M = M+I 
GC TO 80 
I L ( M ) =K 
I U ( M ) = J 
J = L 
F = F+1 
GC TO 30 



GO TC 70 
T) GO TO 20 

T) GO TO 40 
T) GO TO 40 

T) GO TO 40 

T) GO TO 50 
GO TO 30 
J-K) GO TO 60 



70 F=F-L 

1 F ( F .EG. 0) RETURN 
1 = I L ( M) 

J = I U ( M ) 

80 I F ( J- 1 .GE. I I ) GO TO 10 
I F ( I .EG. II i GO TC 5 
1 = 1-1 
90 1=1+1 

I F ( I .EC. J) GO TO 70 
T = A (1 + 1 } 

I F ( A { I ) .LE. Ti GC TO 90 
K = I 

LOO A ( K + 1 ) = A ( K ) 

K — k — 

I F ( T .LT. A ( K ) 5 GO TO 100 
A (K + l ) = T 
GC TO 90 
END 
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